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:

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?

• Any particular reason to use filtfilt instead of filter?
– MBaz
Apr 25, 2016 at 12:34
• There is no particular reason. In fact, I am not sure if filtfilt is even the correct call to use in this case. Apr 25, 2016 at 12:38
• As a general rule, you shouldn't use filtfilt unless you know exactly why you need it.
– MBaz
Apr 25, 2016 at 13:51
• This question is very close to a programming question. I believe a significant problem is identified by the first part of user20703's answer. But the second part of that answer is, I believe, incorrect. keith's answer seems to be the way to go.
– Peter K.
Apr 25, 2016 at 14:09
• @Merin @Merin : your line s = sNew; will zero out your input signal. Why is that there?
– Peter K.
Apr 26, 2016 at 12:46

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

• Thank you very much for this explanation. But do you know what would be the code that can sum all the filter outputs together to form the final output? Apr 25, 2016 at 13:15
• @Merin There's an operator called +… It even works on vectors… Apr 25, 2016 at 13:20
• Thank you so much but I am still getting an error when I try to run this code. The error is in the first line sNew = zeros(s). It says: Error using zeros; Size vector must be a row vector with real elements. Furthermore, is it not better to use the filter function instead of filtfilt? (since we intend our final design to be applied continuously) Apr 25, 2016 at 20:38
• @MarcusMüller I have used this code by keith, and I changed the first line to sNew = zeros(size(s));. So why is the code still not working? I get sNew=0 and the resulting sound is mute. Apr 26, 2016 at 8:16
• @Merin : your line s = sNew; will zero out your input signal. Why is that there?
– Peter K.
Apr 26, 2016 at 12:46

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.

• I made all the changes you've mentioned to my code, I have added my new code to my original post. So, why do I still get sNew = 0? The resulting .wav file is still mute. What is wrong with the code? Apr 26, 2016 at 8:13