Designing a Bandpass Finite Impulse Response (FIR) Filter in Matlab

6 min readNov 20, 2021

In this blog, we will see how to design a non-recursive finite impulse response bandpass filter in Matlab by using the windowing method in conjunction with Kaiser Window.

Introduction

Filter design is an important part of digital signal processing. Filtering is the removal of unwanted frequency components from a signal. The applications of filters are very diverse including audio processing, image processing, telecommunications, etc. Finite impulse response filters are ones’ where
the length of the impulse response is finite. Here we discuss a non-recursive way of realizing a bandpass filter. Bandpass filters pass the frequency components within a certain range of frequencies and attenuate greatly the frequency components lower and higher than the range.

Designing the Filter

Let’s suppose that the following are the filter specifications.

Defining the specifications in Matlab as follows,

%Design of FIR Digital Filter
clc;
close all;
%parameters
Ap = 0.14; %dB maximum passband ripple
Aadesirable = 52; % dB minimum stopband attenuation
wp1 = 500; %rad/s lower passband edge
wp2 = 800; %rad/s upper passband edge
wa1 = 400; % rad/s lower stopband edge
wa2 = 950; %rad/s upper stopband edge
ws = 2400; %rad/s sampling frequency
T = 2*pi/ws;

We can design the filter using the following steps.

Step 1: Selecting the transition width

Selecting the cutoff points

Bt = min(wp1-wa1, wa2-wp2);
wc1 = wp1 - Bt/2;
wc2 = wp2 + Bt/2;

Step 2: Choose δ (represented in code as d) such that the actual passband ripple, Ap is equal to or less than the specified passband ripple⁡ Ãp, and the actual minimum stopband attenuation, A a is equal or greater than
the specified minimum stopband attenuation, Ãa

dp = ((10^(0.05*Ap))-1)/((10^(0.05*Ap))+1);
da = 10^(-0.05*Aadesirable);
d = min(dp, da);

Step 3: With the required δ defined, the actual stopband attenuation Aa can be calculated as

Aa = -20*log10(d); %actual Aa of the filter

Step 4: Choose parameter α (represented in code as alpha) as

%kaiser window
if Aa<=21
   alpha = 0;
elseif Aa<=50
   alpha = (0.5842*((Aa-21)^0.4)) + (0.07886*(Aa-21));
else
   alpha = 0.1102*(Aa - 8.7);
end

Step 5: Choose parameter D as

if Aa<=21
   D = 0.9222;
else
   D = (Aa - 7.95)/14.36;
end

Then select the lowest odd value of N that would satisfy the inequality

N = ceil((ws*D/Bt)+1);
if mod(N,2) == 0
   N=N+1;
end

Form wk (nT) using the following equations:

wk = zeros(N,1);
for n = -(N-1)/2:(N-1)/2
   beta = alpha * (1 - (2*n/(N-1))^2)^0.5;
   numerator = myBessel(beta);
   denominator = myBessel(alpha);
   wk(n+(N-1)/2+1) = numerator/denominator;
end
stem(wk);

A Bessel function of the first kind was written in Matlab as a function for the above purpose. The infinite sum was reduced as follows. It was seen that the summation terms converge to zero. So the terms smaller than 10^(-6) were neglected, thus making it a computable sum.

function [ result ] = myBessel( x )
    k=1;
    result=0;
    term = 10;
    while (term>10^(-6))
        term = (((x/2)^k)/(factorial(k)))^2;
        result = result + term;
        k = k+1;
    end
    
    result = result+1;end

Step 6: Now we have wk (nT). We need h(nT)

w k (nT and h(nT) generated were element-wise multiplied to generate the filter response as follows.

h = zeros(N,1);
h(38) = (2/ws)*(wc2-wc1);
for n = -(N-1)/2:(N-1)/2
   if n==0
      h(n+(N-1)/2+1) = (2/ws)*(wc2-wc1);
   else
      h(n+(N-1)/2+1) = (1/(n*pi)) * (sin(wc2*n*T) - sin(wc1*n*T));
   end
end
figure;
stem(h);

Let’s plot the impulse response of our filter.

fil = h.*wk;
filim =figure;
stem(fil);
title('Impulse Response');
xlabel('n + (N-1)/2');
ylabel('h[n]');
grid on;

Now, we will see how the magnitude response in passband looks like.

[amp, digiFreq] = freqz(fil);analogFreq = digiFreq*ws/(2*pi);
ampdb = 20*log10(abs(amp));fr = figure;
plot(analogFreq, ampdb);
axis([wp1 wp2 -0.05 0.05]);
title('Magnitude Response in Passband');
xlabel('w in rad/s');
ylabel('|H(w)|');
grid on;

Now, we will test the filtering behavior of our filter using an input signal. The input signal given below was generated by the addition of three sin waves.

where ω 1 is the middle frequency of the lower stopband, ω 2 is the middle frequency of the passband, and ω 3 is the middle frequency of the upper stopband

x = zeros(300,1);
w1 = wa1/2;
w2 = (wp2+wp1)/2;
w3 = (ws/2+wa2)/2;for n = 1:300
    x(n) = sin (w1*n*T) + sin (w2*n*T) + sin (w3*n*T);
end

Let’s plot the input signal, expected output signal, and the actual output signal. To find the output fft’s of the input and the filter response were multiplied and ifft was taken.

lenin = length(x);
lenh = length(fil);
lenfft = lenin+lenh-1;IN = fft(x, lenfft);
H = fft(fil, lenfft);
OUT = H.*IN;
out = ifft(OUT, lenfft);threeplots = figure;
subplot(3,1,1);
stem(x);
title('Input Signal');
xlabel('n');
ylabel('x[n]');
grid on;e = zeros(300,1);for n = 1:300
    e(n) = sin (w2*n*T);
end
subplot(3,1,2);
stem(e);
axis([0 300 -2 2]);
title('Expected Output');
xlabel('n');
ylabel('e[n]');
grid on;subplot(3,1,3);
stem(out);
axis([0 300 -2 2]);
title('Output');
xlabel('n');
ylabel('Y[n]');
grid on;

Now, we will see how these three signals behave in the frequency domain.

w = ws*(1-lenfft/2:lenfft/2)/lenfft;
IN1 = abs(fftshift(IN));
OUT1 = abs(fftshift(OUT));
E = fft(e, lenfft);
E1 = abs(fftshift(E));
figure;
subplot(3,1,1);
plot(w, IN1);
title('Frequency Spectrum of Input Signal');
xlabel('w in rad/s');
ylabel('|X(w)|');
grid on;subplot(3,1,2);
plot(w, E1);
title('Frequency Spectrum of Expected Output Signal');
xlabel('w in rad/s');
ylabel('|E(w)|');
grid on;subplot(3,1,3);
plot(w, OUT1);
title('Frequency Spectrum of Output Signal');
xlabel('w in rad/s');
ylabel('|Y(w)|');
grid on;

Discussion

In figure 1 it can be seen that the frequency spectrum of the filter resembles a bandpass filter. Figure 3 shows the ripples in the passband even though they are invisible in figure 1. In figure 5 the first image shows the frequency spectrum of the input signal which is a combination of 3 sinusoids. The second image shows the frequency spectrum of the expected output where only the signal within the passband gets passed to the output. Image 3 shows the actual filtered output which closely resembles the expected output.

Summary and Conclusion

The filter designed for the purpose of this assignment truly acts as a bandpass filter of the stipulated parameters. It is evident from the results shown in figure 5 that the filter blocks the undesired signals and passed the signals within the passband. There is a slight difference between ideally expected
output and the obtained one due to the following reasons.

Passband ripples cause undesirable modifications to the required part of the spectrum.
Non-zero magnitude response in the stopband for the designed filter.
Stopband ripples distort the signals that need attenuation.

You can find the complete code for the filter design at https://github.com/niruhan/fir-filter-design

Acknowledgment: This content was originally submitted as partial fulfillment for the module EN2570: Digital Signal Processing taught by Dr. Chamira Edussooriya in 2016 at the Department of Electronic and Telecommunication Engineering, University of Moratuwa

References

[1] “DESIGN OF NONRECURSIVE (FIR FILTERS)” at http://www.d-
filter.ece.uvic.ca/SupMaterials/Slides/DSP-Ch09-S3,4.pdf.
[2] “ Bessel function ” at https://en.wikipedia.org/wiki/Bessel_function