Implementation of filter bank using fir1

I'm trying to implement a filter bank using fir1 and a kaiser filter.

I have currently implemented my fir1 filter for the lowest frequency band

lpf = fir1(L-1,500/4000,kaiser(L,3));

and now need to shift this filter up to create a total of 10 bands at 250Hz intervals.

How can I shift up the bandpass filter? I would be looking to use some kind of for loop e.g. to apply some function to lpf in order to shift it up by the desired frequency?

Answer 1

[Paul]

Take advantage of the frequency shift property of the DTFT.

For example

Fs = 1000;
f = (-0.5:.001:0.5)*Fs;
L = 50;
lpf = fir1(L-1,500/4000,kaiser(L,3));
h0 = freqz(lpf,1,f,Fs);
h250 = freqz(lpf.*exp(1j*250/Fs*2*pi*(0:(numel(lpf)-1))),1,f,1000);
figure
plot(f,abs([h0(:) h250(:)])),grid
xline(250)
xlabel('Frequency (Hz)');
ylabel('Magnitude')

This approach does result in the shifted filter having complex coefficients, which might not be what you want. If it's not what you want, the problem might have to be reformulated.

[William Greenway]

Thanks very much Paul,

Why do we use 250/Fs rather than just 250? Also I changed Fs to 8000Hz and this no longer seems to work?

[Paul]

The sequence lpf[n] has to be multiplied by exp(1j*wshift*n). So the units of wshift are rad/sample. We have

wshift (rad/sample) = fshift (cycles/sec) / Fs (samples/sec) * (2*pi) (rad/cycle)

It didn't work for Fs = 8000 because I inadvertently hardwired the value of Fs in the call to freqz. Here is the corrected code

Fs = 8000;
f = (-0.5:.001:0.5)*Fs;
L = 50;
lpf = fir1(L-1,500/4000,kaiser(L,3));
h0 = freqz(lpf,1,f,Fs);
h250 = freqz(lpf.*exp(1j*250/Fs*2*pi*(0:(numel(lpf)-1))),1,f,Fs);
figure
plot(f,abs([h0(:) h250(:)])),grid
xline(250)
xlabel('Frequency (Hz)');
ylabel('Magnitude')

[William Greenway]

Thanks very much for your help - it is really helping me to understand this topic! I have managed to get it working!

Final question! I'm finding that when a filter is built, it also has an equivalant negative frequency filter that is reflected in the y axis. When the filter is shifted then the relfected filter can show up. For example see below:

fsamp = 8000;
f = (-0.5:.001:0.5)*fsamp;

% fcuts = [10 30 230  250];
fcuts = [1000 1300 2210 2410];
mags = [0 1 0];
devs = [0.01 0.05 0.01];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fsamp);

n = n + rem(n,2);
lpf = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');

[H0,f0] = freqz(lpf,1,1024,fsamp);
[H250,f250] = freqz(lpf.*exp(1j*2000/fsamp*2*pi*(0:(numel(lpf)-1))),1,1024,fsamp);
plot(abs([f0 f250]),abs([H0 H250]))

As can be seen by the orange line, when the filter is shifted, the negative reflection of the filter is brought up into the positive frequency domain. How can I get rid of the negative frequency components of lpf?

[Paul]

The short anwer to the the question is: You can't, at least not without further manipulation to achieve whatever the goal is.

The reason is that the frequency response of a discrete time filter is a periodic function defined over the entire real frequency line, with period 2*pi (rad) or (2*pi)/Ts (rad/sec) or Fs (Hz) depending on which frequency you prefer. Typically, the frequecy response is only plotted over a half period, but we can see the full picture by plotting the magnitude response of lpf over three periods (we'll only plot the magnitude here, but to get the full picture you should plot the phase as well). Note also that lpf is now a band pass filter, rather than a low pass filter as in the original code.

fsamp = 8000;
f = (-0.5:.001:0.5)*fsamp;
% fcuts = [10 30 230  250];
fcuts = [1000 1300 2210 2410];
mags = [0 1 0];
devs = [0.01 0.05 0.01];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fsamp);
n = n + rem(n,2);
lpf = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
[H0,f0] = freqz(lpf,1,3*f(:),fsamp);
figure
plot(f0,abs(H0))

So when the frequency reponse of the filter is shifted right by 2000 Hz, as you said a portion of the first band at negative frequency will be shifted up into low positive frequency, and the first band at positive frequency will be shifted up so a portion of the shifted band is cut off at the Nyquist frequency (Fs/2).

figure;
plot(f0+2000,abs(H0),'r')
xlim([-1.5 1.5]*1e4) % make the x axis limits the same as the previous plot
xline(0)
xline(Fsamp/2)

If the frequency response of a filter is shifted, its periodicity will have the observed effect.

[William Greenway]

Sure that makes sense! Thanks very much for your answer :)

I'm trying to design a bank of 10 filters shifted by 250Hz each filter. Is there any way of creating a bank of 10 filters without the periodicity having the effect that you mentioned (i.e. getting rid of the other periods than the first half period)?

Also why is lpf now a bandpass filter when before it was a low pass filter?

[Paul]

As to the second question, it looks like the filter design parameters changed. As shown, lpf from the original code was low pass and from the new code is bandpass.

Fs = 8000;
f = (-0.5:.001:0.5)*Fs;
L = 50;
lpf = fir1(L-1,500/4000,kaiser(L,3));  % lpf in the original question

fcuts = [1000 1300 2210 2410];
mags = [0 1 0];
devs = [0.01 0.05 0.01];
fsamp = Fs;
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fsamp);
n = n + rem(n,2);
lpfnew = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale'); % lpf from comment above

[h0,f0] = freqz(lpf,1,1024,Fs);
[hnew,fnew] = freqz(lpfnew,1,1024,Fs);
figure
plot([f0 fnew],abs([h0 hnew])),grid
legend('lpf','lpfnew')

As to the first question, maybe you could achieve your goal by designing 10 different filters and adding them together (rather than shifted copies of a single filter) to get the desired, overall frequency response.