Lab Handout 09: FIR Filter Design Using MATLAB
Course Code: EEE 420Course Name: Digital Signal Processing
Electrical & Electronic Engineering , Eastern Mediterranean University
Instructor: Asst. Prof. Dr. Erhan A. İnce
E-mail:
Ext: 2778 / Lab Assistant: Mr. Burçin Özmen
E-mail:
Ext: 2771
LINEAR-PHASE FIR FILTER
Shapes of impulse and frequency responses and locations of system function zeros of linear-phase FIR filters will be discussed.
Let , be the impulse response of length (or duration) M. Then the system function is
which has poles at the origin z = 0 (trivial poles) and zeros located anywhere in the z-plane. The frequency response function is
Type-1 Linear-Phase FIR Filters
Symmetrical impulse response and M odd. In this case , is an integer, and , . Then we can show that
where sequence is obtained from as
Since then, we have
where is an amplitude response function and not a magnitude response function.
function [Hr,w,a,L]=Hr_Type1(h);
% Computes Amplitude response Hr(w)of a Type-1 LP FIR filter
% ------
% [Hr,w,a,L]=Hr_Type1(h)
% Hr = Amplitude Response
% w = 500 frequencies between [0,pi] over which Hr is computed
% a = Type-1 LP filter coefficients
% L = Order of Hr
% h = Type-1 LP filter impulse response
%
M = length(h);
L= (M-1)/2;
a = [h(L+1) 2*h(L:-1:1)]; % 1x(L+1) row vector
n = [0:1:L] % (L+1)x1 column vector
w = [0:1:500]'*pi/500;
Hr = cos(w*n)*a';
Type-2 Linear-Phase FIR Filters
Symmetrical impulse response and M even. In this case , is not an integer, and , . Then we can show that
where
Hence,
where is an amplitude response function. Note that we get
regardless of or . Hence we cannot use this type (i.e. symmetrical and M even.) for highpass or bandstop filters.
Type-3 Linear-Phase FIR Filters
Antisymmetric impulse response and M odd. In this case , is an integer, and , , and . Then we can show that
where
Hence,
Note that at and we have , regardless of or . Furthermore, , which means that is purely imaginary. Hence this type of filter is not suitable for designing a lowpass or a highpass filter. However, this behavior is suitable for approximating ideal Hilbert transformers and differentiators.
function [Hr,w,c,L]=Hr_Type3(h);
% Computes Amplitude response Hr(w) of a Type-3 LP FIR filter
% ------
% [Hr,w,c,L]=Hr_Type3(h)
% Hr = Amplitude Response
% w = frequencies between [0,pi] over which Hr is computed
% b = Type-3 LP filter coefficients
% L = Order of Hr
% h = Type-3 LP filter impulse response
%
M = length(h);
L = (M-1)/2;
c = [2*h(L+1:-1:1)];
n = [0:1:L];
w = [0:1:500]'*pi/500;
Hr = sin(w*n)*c';
Type-4 Linear-Phase FIR Filters
Antisymmetric impulse response and M even. This case is similar to Type-2. We have
where
and
Note that at and we have and . Hence, this type is also suitable for designing digital Hilbert transformers and differentiators.
Assignments
- MATLAB functions for Type-1 and Type-3 is given. Write a MATLAB function to compute the amplitude response and the location of the zeros of for Type-2 and Type-4.
- Let
Determine the Type of the Linear-Phase FIR filters. Find the amplitude response and locations of the zeros of . Also plot the impulse response, amplitude response, coefficients, and the zero locations.
WINDOW DESIGN TECHNIQUES
Summary of commonly used window function characteristics.
WindowName / Transition Approximate / Width Exact Values / Min. Stopband Attenuation
Rectangular / / / 21 dB
Bartlett / / / 25 dB
Hanning / / / 44 dB
Hamming / / / 53 dB
Blackman / / / 74 dB
MATLAB provides several routines to implement window functions discussed above table.
· w = boxcar(n) returns the n-point rectangular window in array w.
· w = triang(n) returns the n-point Bartlett (Triangular) window in array w.
· w = hanning(n) returns the n-point symmetric Hanning window in a column vector in array w .
· w = hamming(n) returns the n-point symmetric Hamming window in a column vector in array w.
· w = blackman(n) returns the n-point symmetric Blackman window in a column vector in array w.
· w = kaiser(n,beta) returns the beta-valued n-point Kaiser window in array w.
Using these routines, we can use MATLAB to design FIR filters based on the window technique, which also requires an ideal lowpass impulse response as shown below.
function hd=ideal_lp(wc,M);
% Ideal LowPass filter computation
% ------
% [hd] = ideal_lp(wc,M);
% hd = ideal impulse response between 0 to M-1
% wc = cutoff frequency in radians
% M = length of the ideal filter
%
alpha = (M-1)/2;
n = [0:1:(M-1)];
m = n - alpha +eps; % add smallest number to avoi divided by zero
hd = sin(wc*m)./(pi*m);
To display the frequency domain plots of digital filters, MATLAB provides the freqz routine. Using this routing, we can developed a modified version, called freqz_m, which returns the magnitude response in absolute as well as dB scale, the phase response, and the group delay response as shown below.
function [db,mag,pha,grd,w] = freqz_m(b,a)
% Modified version of freqz subroutine
% ------
% [db,mag,pha,grd,w] = freqz_m(b,a)
% db = relative magnitude in dB computed over 0 to pi radians
% mag = absolute magnitude computed over 0 to pi radians
% pha = Phase response in radians over 0 to pi radians
% grd = Group delay over 0 to pi radians
% w = 501 frequency samples between 0 to pi radians
% b = numerator polynomial of H(z) (for FIR: b=h)
% a = denominator polynomial of H(z) (for FIR: a=[1])
%
[H,w] = freqz(b,a,1000,'whole') ;
H = (H(1:1:501))'; w = (w(1:1:501))';
mag = abs(H);
db = 20*log10((mag+eps)/max(mag));
pha = angle(H);
grd = grpdelay(b,a,w);
Assignment 3
- Design a digital FIR lowpass filter with the following specifications:
Choose an appropriate window function from above Table. Determine the impulse response and provide a plot of the frequency response of the designed filter.
Example 1: Let us design the following digital bandpass filter.
These quantities are shown in the following figure.
% Solution of Example 1
ws1 = 0.2*pi; wp1 = 0.35*pi;
wp2 = 0.65*pi; ws2 = 0.8*pi;
As = 60;
tr_width = min((wp1-ws1),(ws2-wp2));
M = ceil(11*pi/tr_width) +1;
n = [0:1:M-1];
wc1 = (ws1+wp1)/2; wc2 = (wp2+ws2)/2;
hd = ideal_lp(wc2,M) - ideal_lp(wc1,M);
w_bla = (blackman(M))';
h = hd.*w_bla;
[db,mag,pha,grd,w] = freqz_m(h,[1]);
delta_w =2*pi/1000;
Rp = -min(db(wp1/delta_w+1:1:wp2/delta_w)) % Actual Passband Ripple
As = -round(max(db(ws2/delta_w+1:1:501))) % Min Stopband Attenuation
% Plots
subplot(2,2,1); stem(n,hd); title('Ideal Impulse Response')
axis([0 M-1 -0.4 0.5]); xlabel('n'); ylabel('hd(n)')
subplot(2,2,2); stem(n,w_bla); title('Blackman Window')
axis([0 M-1 0 1.1]); xlabel('n'); ylabel('w(n)')
subplot(2,2,3); stem(n,h); title('Actual Impulse Response')
axis([0 M-1 -0.4 0.5]); xlabel('n'); ylabel('h(n)')
subplot(2,2,4); plot(w/pi,db); axis([0 1 -150 10]);
title('Magnitude Response in db');grid;
xlabel('frequency in pi units'); ylabel('Decibels')
Assignment 4
- The frequency response of an ideal bandstop filter is given by
using a Kaiser window, design a bandstop filter of length 45 with stopband attenuation of 60 dB.
Example 2: The frequency response of an ideal digital differentiator is given by
using a Hamming window of length 21, design a digital FIR differentiator. Plot the time- and the frequency-domain responses.
Soln: Note that M is an even number, then is an integer and will be zero for all n. hence M must be an odd number, and this will be a Type-3 linear phase FIR filter.
% Solution of Example 2
M = 21; alpha = (M-1)/2;
n = 0:1:M-1;
hd = (cos(pi*(n-alpha)))./(n-alpha); hd(alpha+1) = 0;
w_ham = (hamming(M))';
h = hd.*w_ham;
[Hr,w,P,L]=Hr_Type3(h);
% Plots
% subplots(1,1,1);
subplot(2,2,1); stem(n,hd); title('Ideal Impulse Response')
axis([-1 M -1.2 1.2]); xlabel('n'); ylabel('hd(n)')
subplot(2,2,2); stem(n,w_ham); title('Hamming Window')
axis([-1 M 0 1.2]); xlabel('n'); ylabel('w(n)')
subplot(2,2,3); stem(n,h); title('Actual Impulse Response')
axis([-1 M -1.2 1.2]); xlabel('n'); ylabel('h(n)')
subplot(2,2,4); plot(w/pi,Hr/pi); title('Amplitude Response');grid;
xlabel('frequency in pi units'); ylabel('slope in pi units');
axis([0 1 0 1]);