E E 2 7 5 Lab June 30, 2006

FIR Filters in Matlab

Lab 5. FIR Filter Design in Matlab

Digital filters with finite-duration impulse reponse (all-zero, or FIR filters) have both advan-
tages and disadvantages when compared to infinite-duration impulse response (IIR) filters.

FIR filters have the following primary advantages:

• They can have exactly linear phase.
• They are always stable, even when quantized.
• The design methods are generally linear.
• They can be realized efficiently in hardware.
• The filter startup transients have finite duration.

The primary disadvantage of FIR filters is that they often require a much higher filter order
than IIR filters to achieve a given level of performance. Correspondingly, the delay of these
filters is often much greater than for an equal performance IIR filter.

FIR Methods

The Signal Processing Toolbox supports a variety of methods for the design of FIR filters.

F ilterM ethod Description F ilterF unctions

Windowing Apply window to truncated inverse fir1, fir2, kaiserord
Fourier transform of desired filter
Multiband with Equiripple or least squares approach firls, firpm, firpmord
Transition Bands over frequency subbands
Constrained Least Minimize squared integral error over fircls, fircls1
Squares entire frequency range subject to maximum
error constraints
Arbitrary Response Arbitrary responses, including nonlinear cfirpm
phase and complex filters
Raised Cosine Lowpass response with smooth, firrcos
sinusoidal transition

2006GM
c
Impulse Response Revisited

FIR filters are described by difference equations of the form

M
X
y(n) = bm x(n − m)
m=0

where the filtered signal y(n) is just a linear combination of current and previous values
of the input signal x(n). The coefficients b are the numerator coefficients of the transfer
function. The denominator of the transfer function will always be a = 1. The order of the
filter is n = length(b) − 1.

If the input signal is the unit impulse x = [1 0 0 0 ...], then the corresponding impulse response
y(n) = h(n) is identical to b(n):

h(0) = b0 x(0) = b0
h(1) = b0 x(1) + b1 x(0) = b1
h(2) = b0 x(2) + b1 x(1) + b2 x(0) = b2 ...etc.

The FIR filter coefficients give the impulse response.

Try this:
b = [-1 0 2 -3 1]; % no need to specify the a coefficients
stem(b)
figure, impz(b,1)

Why are the x-axis scales different?

Linear Phase Filters

A filter whose impulse response is symmetric about its midpoint is called a (generalized)
linear phase filter. For such filters,
• The DFT of the impulse response will be either purely real or purely imaginary.
• The magnitude of the DFT is scaled by the filter’s magnitude response (there is no ampli-
tude distortion).
• The phase shift φ of a filtered signal will vary linearly with frequency ω (pure time delay
with no phase distortion).
• The phase delay −φ(ω)/ω and group delay −dφ(ω)/d(ω) will be equal and constant. For
an order n linear phase FIR filter, the phase delay and group delay is n/2.
The absence of either amplitude distortion or phase distortion preserves the waveform of
signals in the passband.

Except for cfirpm, all the FIR filter design functions in the Signal Processing
Toolbox design linear phase filters only.

Try:
a = 1;
b = fir1(5,0.5);
fvtool(b,a)

Look at the phase delay and the group delay.

FIR Filter Types

The symmetric impulse response of a linear phase filter can have an odd or an even number
of points,and can have an odd or even symmetry about the midpoint, leading to four filter
types:

• Type I: Odd length, even symmetry

• Type II: Even length, even symmetry
• Type III: Odd length, odd symmetry
• Type IV: Even length, odd symmetry

Depending on the filter type, certain restrictions apply:

F ilter F ilter Symmetry ResponseH(0) ResponseH(1)

T ype Order (N yquist)
Type I Even b(k) = b(n + 2 − k), k = 1, ..., n + 1 No restriction No restriction
Type II Odd b(k) = b(n + 2 − k), k = 1, ..., n + 1 No restriction H(1) = 0
Type III Even b(k) = −b(n + 2 − k), k = 1, ..., n + 1 H(0) = 0 H(1) = 0
Type IV Odd b(k) = −b(n + 2 − k), k = 1, ..., n + 1 H(0) = 0 No restriction

The functions fir1, fir2, firls, firpm, fircls, fircls1, and firrcos all design type I and II linear
phase FIR filters by default. Both firls and firpm design type III and IV linear phase FIR
filters given a ’hilbert’ or ’differentiator’ flag. The function cfirpm can design any type of
linear or nonlinear phase filter.
Because the frequency response of a type II filter is zero at the Nyquist frequency (“high”
frequency), fir1 does not design type II highpass and bandstop filters. For odd-valued n in
these cases, fir1 adds 1 to the order and returns a type I filter.

Window-Based Design

Windowing is a common design method for FIR filters. In this method,

• The ideal frequency response H(f ) is sampled.
• The corresponding ideal impulse response h(n) is determined by the inverse Fourier trans-
form. In general, this response cannot be implemented in a digital filter because it is infinite
and noncausal.
• h(n) is symmetrically truncated (multiplied by a finite, symmetric window) to create a
linear phase finite impulse response.

The approximation to the ideal filter is “best” in a mean square sense, compared to other
approximations of the same length, by Parseval’s theorem. However, the abrupt truncation
leads to overshoot (Gibb’s phenomenon) and ripples in the spectrum. The undesirable effects
of truncation are reduced or eliminated by the use of tapered windows.

Windowing does not explicitly impose amplitude response constraints, such as passband
ripple or stopband attenuation. It must be used iteratively to produce designs that meet
such specifications.

Try:
edit windemo
windemo(20)
windemo(50)
windemo(100)

What do the arguments (20, 50, 100) specify?

How does the filter change with the argument?

Windowing Functions

The Signal Processing Toolbox supports a variety of windows commonly used in FIR filter
design. Typing
help window
provides a list of available functions:
bartlett Bartlett window
barthannwin Modified Bartlett-Hanning window
blackman Blackman window
blackmanharris Minimum 4-term Blackman-Harris window
bohmanwin Bohman window
chebwin Chebyshev window
flattopwin Flat Top window
gausswin Gaussian window
hamming Hamming window
hann Hann window
kaiser Kaiser window
nuttallwin Nuttall defined minimum 4-term Blackman-Harris window
parzenwin Parzen (de la Valle-Poussin) window
rectwin Rectangular window
triang Triangular window
tukeywin Tukey window

Individual functions take inputs for a window length n and window parameters and return
the window w in a column vector of length n. (Note: Use w’ for array products with row vec-
tor impulse responses.) The window function serves as a gateway to the individual functions.

w = gausswin(64,alpha)
and
w = window(@gausswin,64,alpha)
both return a Gaussian window of length 64 with standard deviation equal to 1/alpha.

Try:
help window
n=15;
w=gausswin(n,3);
stem(-7:7,w)
n=linspace(0,6*pi,100);
y=sinc(n);
figure, stem(n,y)
W=gausswin(2*length(y),3);
w=W(length(y)+1:end);
wy=w’.*y;
figure, stem(n,wy,’m’)
Explain the y and wy plots.

Windowing and Spectra

When a signal is truncated, high-frequency components are introduced that are visible in
the DFT. By windowing the truncated signal in the time domain, endpoints are assigned a
reduced weight. The effect on the DFT is to reduce the height of the side lobes, but increase
the width of the main lobe.

Try:
edit windft
windft

Tuncated signal and DFT:

Truncated Signal
1

0.5

−0.5

−1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

DFT of Truncated Signal

50

40
Magnitude

30

20

10

0
0 1 2 3 4 5 6 7 8 9 10

Tuncated, windowed signal and DFT:

 Truncated, Windowed Signal
1

0.5

−0.5

−1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

DFT of Truncated, Windowed Signal

25

20
Magnitude

15

10

0
0 1 2 3 4 5 6 7 8 9 10

Window Visualization Tool

The transition bandwidth of a window-based FIR filter is determined by the width of the
main lobe of the DFT of the window function, adjustable by changing the filter order.

Passband and stopband ripples are determined by the magnitude of the side lobe of the DFT
of the window function, and are usually not adjustable by changing the filter order.

The actual approximation error is scaled by the amount of the passband magnitude response.

Ideally, the spectrum of a window should approximate an impulse. The main lobe should be
as narrow as possible and the side lobes should contain as little energy as possible.

The Window Visualization Tool (WVTool) allows you to investigate the tradeoffs among
different windows and filter orders.

wvtool(windowname(n))
opens WVTool with time and frequency domain plots of the n-length window specified in
windowname, which can be any window in the Signal Processing Toolbox. Several windows
can be given as input arguments for comparative display.
Try:
wvtool(hamming(32),kaiser(32,2.5),flattopwin(32))
wvtool(kaiser(32,1),kaiser(32,5),kaiser(32,10))

Window Design and Analysis Tool

The Window Design and Analysis Tool (WinTool) is used in conjunction with the Window
Visualization Tool.

• Use WVTool for displaying and comparing existing windows created in the Matlab workspace.
• Use WinTool to interactively design windows with certain specifications and export them
to the Matlab workspace.

Most window types satisfy some optimality criterion. Some windows are combinations of
simpler windows. For example, the Hann window is the sum of a rectangular and a cosine
window, and the Bartlett window is the convolution of two rectangular windows. Other
windows emphasize certain desirable features. The Hann window improves high-frequency
decay (at the expense of larger peaks in the side lobes). The Hamming window minimizes
side lobe peaks (at the expense of slower high-frequency decay). The Kaiser window has
a parameter that can be tuned to control side lobe levels. Other windows are based on
simple mathematical formulas for easy application. The Hann window is easy to use as a
convolution in the frequency domain.

An optimal time-limited window maximizes energy in its spectrum over a given frequency
band. In the discrete domain, the Kaiser window gives the best approximation to such an
optimal window.

wintool
opens WinTool with a default 64-point Hamming window. Try it - experiment with different
window designs and export them to the workspace.

Comment on the Chebyshev compared to the Blackman-Harris window.

Can you think of an advantage one would have over the other?
Example: Lowpass Filter

Consider the ideal, or “brick wall”, digital low-pass filter with a cutoff frequency of ω0 rad/s.
This filter has magnitude 1 at all frequencies less than ω0 , and magnitude 0 at frequencies
between ω0 and π. Its impulse response sequence h(n) is
1 Zπ 1 Z ω0 jωn ω0 ω0
 
jωn
H(ω)e dω = e dω = sinc n
2π −π 2π −ω0 π π

This filter is not implementable since its impulse response is infinite and noncausal. To
create a finite-duration impulse response, turncate it by applying a window. Retain the
central section of impulse response in the turncation to obtain a linear phase FIR filter. For
example, a length 51 filter with a lowpass cutoff frequency ω0 of 0.4π rad/s is
b=0.4*sinc(0.4*(-25:25));

The window applied here is a simple rectangular window. By Parseval’s theorem, this is
the length 51 filter that best approximates the ideal lowpass filter, in the integrated least
squares sense. To display the filter’s frequency response in FVTool, type
fvtool(b,1)

Ringing and ripples occur in the response, especially near the band edge. This “Gibb’s ef-
fect” does not vanish as the filter length increases, but a nonrectangular window reduces its
magnitude. Multiplication by a window in the time domain causes a convolution or smooth-
ing in the frequency domain. Apply a length 51 Hamming window to the filter and display
the result using FVTool:
bw=b.*hamming(51)’;
fvtool(b,1)

Using a Hamming window greatly reduces the ringing. This is at the expense of transition
width (the windowed version takes longer to ramp from passband to stopband) and optimal-
ity (the windowed version does not minimize the integrated least squared error).

Try:
b=0.4*sinc(0.4*(-25:25));
fvtool(b,1)
bw=b.*hamming(51)’;
fvtool(bw,1)
Right-click the y-axis label in FVTool and choose Magnitude squared on both plots.
Where in the plot is the ringing reduced by the window?

Standard Band FIR Design

The Signal Processing Toolbox functions fir1 and fir2 are both based on the windowing
method. Given a filter order and a description of an ideal filter, these functions return a
windowed inverse Fourier transform of the ideal filter. Both use Hamming windows by de-
fault, but they accept any windowing function.

fir1 resembles the IIR filter design functions in that it is formulated to design filters in stan-
dard band configurations: lowpass, bandpass, highpass, and bandstop. The commands
n=50;
Wn=0.4;
b=fir1(n,Wn);
create a row vector b containing the coefficients of the order n Hamming-windowed filter.
This is a lowpass, linear phase FIR filter with cutoff frequency W n. W n is a number between
0 and 1, where 1 corresponds to the Nyquist frequency, half the sampling frequency.

For a highpass filter, simply append the string ’high’ to the function’s parameter list. For
a bandpass or bandstop filter, specify Wn as a two-element vector containing the passband
edge frequencies; append the string ’stop’ for the bandstop configuration.

b=fir1(n,Wn,window);
uses the window specified in column vector window for the design. The vector window must
be n + 1 elements long. If you do not specify a window, fir1 applies a Hamming window.

The kaiserord function estimates the filter order, cutoff frequency, and Kaiser window β pa-
rameter needed to meet a given set of specifications. Given a vector of frequency band edges,
a vector of magnitude, and a maximum allowable ripple, kaiserord returns appropriate input
parameters for the fir1 function.

Try:
edit fir1demo
fir1demo

Exercise
 1. Design a windowed FIR bandstop filter to remove the 300 Hz component from the three-
tone signal with noise y from Lab IIR Filters in Matlab. Use a sampling frequency of
8192 Hz.

2. Plot the response.

3. Filter the signal y with the designed filter.

4. Compare signals and spectra before and after filtering.

Arbitrary Response FIR Filters

The fir2 function also designs windowed FIR filters, but with an arbitrarily shaped piecewise
linear frequency response. (The IIR counterpart of this function is yulewalk).

b=fir2(n,f,m);
returns row vector b containing the n + 1 coefficients of an order n FIR filter. The frequency-
magnitude characteristics of this filter match those given by vectors f and m.

b=fir2(n,f,m,window);
uses the window specified in column vector window for the design. The vector window must
be n + 1 elements long. If you do not specify a window, fir2 applies a Hamming window.

The function cfirpm is used to design complex and nonlinear-phase equiripple FIR filters. It
allows arbitrary frequency-domain constraints.

Try:
edit directstop2
directstop2(10)
directstop2(100)
directstop2(500)
edit nlinphase
nlinphase
Multiband Filters

The function firls designs linear-phase FIR filters that minimize the weighted, integrated
squared error between an ideal piecewise linear function and the magnitude response of the
filter over a set of desired frequency bands.

b=firls(n,f,a)
returns row vector b containing the n+1 coefficients of the order n FIR filter whose frequency-
amplitude characteristics approximately match those given by vectors f and a.

The function firls allows you to introduce constraints by defining upper and lower bounds
for the frequency response in each band.

The function fircls1 is used specifically to design lowpass and highpass linear phase FIR
filters using constrained least squares.

Try:
edit firlsdemo
firlsdemo

edit firclsdemo
firclsdemo

Raised Cosine Filters

The sinc function, which is the impulse response of an ideal lowpass filter, forms the basis
for several other interpolating functions of the form
n
 
h(n) = f (n)sinc , where f (0) = 1
ns
One commonly-used form is the raised cosine function:

cos(πR nns ) n
 
hrc (n) = n 2 sinc , 0≤R≤1
1 − (2R ns ) ns

R is called the rolloff factor. Like the sinc function, the raised cosine function is 1 at n = 0
and 0 at all other sampling instances n = ns . In contrast to the sinc function, the raised
cosine has faster decaying oscillations on either side of the origin for R > 0. This results in
improved reconstruction if samples are not acquired at exactly the sampling instants (i.e., if
there is jitter). It also uses fewer past and future values in the reconstruction, as compared
to the sinc function. The shape of the function’s spectrum is the “raised cosine”.

The ideal raised cosine lowpass filter frequency response consists of unity gain at low fre-
quencies, a raised cosine function in the middle, and total attenuation at high frequencies.
The width of the transition band is determined by the rolloff factor.

b=firrcos(n,F0,df,fs)
b=firrcos(n,F0,df,fs,’bandwidth’)
are equivalent, and return an order n lowpass linear-phase FIR filter with a raised cosine
transition band. The cutoff frequency is F 0, the transition bandwidth df , and sampling
frequency is f s, all in hertz. df must be small enough so that F 0 ± df /2 is between 0 and
f s/2. b is normalized so that the nominal passband gain is always equal to 1.

b=firrcos(n,F0,r,fs,’rolloff’)
interprets the third argument, r, as the rolloff factor instead of the transition bandwidth,
df . r must be in the range [0, 1].

Try:
b = firrcos(20,0.25,0.25,2);
freqz(b)

Describe what you see.

Frequency Domain Filtering

Often, a long (perhaps continuous) stream of data must be processed by a system with only
a finite length buffer for storage. The data must be processed in pieces, and the processed
result constructed from the processed pieces.

In the overlap-add method, an input signal x(n) is partitioned into equal length data blocks.
The filter coefficients (impulse response) and each block of data are transformed to the fre-
quency domain using the FFT, where they can be efficiently convolved using multiplication.
The partial convolutions of the signal are returned to the time domain with the IFFT, where
they are shifted and summed using superposition.
fftfilt implements the overlap-add method for FIR filters.

y=fftfilt(b,x)
uses an FFT length of nfft = 2nextpow2(n) and a data block length of nfft-length (b)+1 (ensures
circular convolution).

firfilt incurs an “offline” startup cost when converting the coefficients b to the frequency
domain. After that, the number of multiplications fftfilt performs relative to filter (which
implements the filter in direct form) is ≈ log2 (L)/N , where L is the block length and N is the
filter length. (Multiplications are a good measure of performance, since they are typically
expensive on hardware.) The net result is that for filters of high order, fftfilt outperforms
filter.

Try:
x=[1 2 3 4 5 6];
h=[1 1 1];
y=conv(h,x)
x1=[1 2 3];
x2=[4 5 6];
y1=conv(h,x1)
y2=conv(h,x2)
Y1=[y1, zeros(1,3)];
Y2=[zeros(1,3),y2];
Y=Y1+Y2

Describe what is happening in this code.

The following script takes a few moments to run.

edit filttimes
filttimes

Above what order can you say that fftfilt is always faster?

(The material in this lab was put together by Paul Beliveau and handout derives principally from
the MathWorks training document “MATLAB for Signal Processing”, 2006.)

2006GM
c