Filter Design in Matlab

Question

I want to design a filter to remove all frequency components from a .wav file except those within $\pm 25\,\text{Hz}$ of $523\,\text{Hz}$ as well as its harmonics (up to the Nyquist frequency).

Attempt

I tried to suppress frequencies below 523-25 Hz using a high-pass filter, and suppress the ones higher than 523+25 Hz using a low pass filter. Here's my code:

[s, Fs] = wavread('x.wav');

for N=2:1:(Fs/(2*(523+25)))
    % High-pass filter
    FsNorm =  (523-25).*N / (Fs/2);
    [b,a] = butter(10, FsNorm, 'high');
    sHigh = filtfilt(b, a, s);

    % Low pass
    FsNorm = (523+25).*N / (Fs/2);
    [b,a] = butter(10, FsNorm, 'low'); 
    sLow = filtfilt(b, a, s);
end

wavwrite(sLow, Fs, 'x_filtered');

The code doesn't work, the resulting file is mute. My intention is to apply each filter to the output of the previous filter, but I don't know how to implement this in code.

I used $(523+25) N_{max} < \frac{F_s}{2}$ to include all the harmonics up to the Nyquist frequency.

So what do I need to do to get this filter to work?

Edit:

Here's my new code:

[s, Fs] = wavread('x.wav');

sNew = zeros(size(s));

for N=2:1:(Fs/(2*(523+25)))
    % High-pass filter
    FsNorm =  (523-25).*N / (Fs/2);
    [b,a] = butter(10, FsNorm, 'high');
    sHigh = filter(b, a, s);

    % Low pass
    FsNorm = (523+25).*N / (Fs/2);
    [b,a] = butter(10, FsNorm, 'low'); 
    sLow = filter(b, a, s);

    Now sum the results
    sNew = sHigh + sLow;
end

sNew

wavplay(sNew,Fs);
wavwrite(sNew, Fs, 'x_filtered');

I still get sNew=0 and the sound is still mute. What's wrong?

Answer 1

You just need to make a simple modification:

Send the original signal through each filter and sum the output of each filter together to get your final result.

At the moment you are passing the filtered data on to each subsequent filter. If you think carefully about that, the output will be silent i.e. you filter out everything but 523 Hz +/- 25 which leaves you just frequencies around 523 Hz, then you take that and filter out everything but 1046 Hz +/- 25 which will strip out the frequencies around 523 Hz and then you have silence after just two iterations of your loop.

---

s = s(:).';
sNew = zeros(1,length(s)); % The eventual output

for N=2:1:(Fs/(2*(523+25)))
    % High-pass filter
    FsNorm =  (523-25).*N / (Fs/2);
    [b,a] = butter(10, FsNorm, 'high');
    sHigh = filter(b, a, s);

    % Low pass
    FsNorm = (523+25).*N / (Fs/2);
    [b,a] = butter(10, FsNorm, 'low'); 
    sLow = filter(b, a, sHigh); % Note we are filtering sHigh, not 's'

    % Now sum the results
    sNew = sNew + sLow;
end

Answer 2

You code currently produces a signal sLow that is s through a single low-pass filter (using the last value of N), because you keep applying all the filtering to the original signal s.

So you probably need to substitute s with sNew in the loop. Before the loop starts do s = sNew. At the end of each loop do sNew = sLow + sHigh. Then sNew is your final result.