Print Preview C:\TEMP\Apdf_2541_3068\home\AppData\Local\PTC\Arbortext\Editor\.aptcache\ae3849j0/tf38450r Signal Processing Toolbox User's Guide

User Manual:

Open the PDF directly: View PDF PDF.
Page Count: 527 [warning: Documents this large are best viewed by clicking the View PDF Link!]

Signal Processing Toolbox™
User’s Guide
R2012b
How to Contact MathWorks
www.mathworks.com Web
comp.soft-sys.matlab Newsgroup
www.mathworks.com/contact_TS.html Technical Support
suggest@mathworks.com Product enhancement suggestions
bugs@mathworks.com Bug reports
doc@mathworks.com Documentation error reports
service@mathworks.com Order status, license renewals, passcodes
info@mathworks.com Sales, pricing, and general information
508-647-7000 (Phone)
508-647-7001 (Fax)
The MathWorks, Inc.
3 Apple Hill Drive
Natick, MA 01760-2098
For contact information about worldwide offices, see the MathWorks Web site.
Signal Processing Toolbox™ User’s Guide
© COPYRIGHT 1988–2012 by The MathWorks, Inc.
The software described in this document is furnished under a license agreement. The software may be used
or copied only under the terms of the license agreement. No part of this manual may be photocopied or
reproduced in any form without prior written consent from The MathWorks, Inc.
FEDERAL ACQUISITION: This provision applies to all acquisitions of the Program and Documentation
by, for, or through the federal government of the United States. By accepting delivery of the Program
or Documentation, the government hereby agrees that this software or documentation qualifies as
commercial computer software or commercial computer software documentation as such terms are used
or defined in FAR 12.212, DFARS Part 227.72, and DFARS 252.227-7014. Accordingly, the terms and
conditions of this Agreement and only those rights specified in this Agreement, shall pertain to and govern
theuse,modification,reproduction,release,performance,display,anddisclosureoftheProgramand
Documentation by the federal government (or other entity acquiring for or through the federal government)
and shall supersede any conflicting contractual terms or conditions. If this License fails to meet the
government’s needs or is inconsistent in any respect with federal procurement law, the government agrees
to return the Program and Documentation, unused, to The MathWorks, Inc.
Trademarks
MATLAB and Simulink are registered trademarks of The MathWorks, Inc. See
www.mathworks.com/trademarks for a list of additional trademarks. Other product or brand
names may be trademarks or registered trademarks of their respective holders.
Patents
MathWorks products are protected by one or more U.S. patents. Please see
www.mathworks.com/patents for more information.
Revision History
1988 First printing New
November 1997 Second printing Revised
January 1998 Third printing Revised
September 2000 Fourth printing Revised for Version 5.0 (Release 12)
July 2002 Fifth printing Revised for Version 6.0 (Release 13)
December 2002 Online only Revised for Version 6.1 (Release 13+)
June 2004 Online only Revised for Version 6.2 (Release 14)
October 2004 Online only Revised for Version 6.2.1 (Release 14SP1)
March 2005 Online only Revised for Version 6.2.1 (Release 14SP2)
September 2005 Online only Revised for Version 6.4 (Release 14SP3)
March 2006 Sixth printing Revised for Version 6.5 (Release 2006a)
September 2006 Online only Revised for Version 6.6 (Release 2006b)
March 2007 Online only Revised for Version 6.7 (Release 2007a)
September 2007 Online only Revised for Version 6.8 (Release 2007b)
March 2008 Online only Revised for Version 6.9 (Release 2008a)
October 2008 Online only Revised for Version 6.10 (Release 2008b)
March 2009 Online only Revised for Version 6.11 (Release 2009a)
September 2009 Online only Revised for Version 6.12 (Release 2009b)
March 2010 Online only Revised for Version 6.13 (Release 2010a)
September 2010 Online only Revised for Version 6.14 (Release 2010b)
April 2011 Online only Revised for Version 6.15 (Release 2011a)
September 2011 Online only Revised for Version 6.16 (Release 2011b)
March 2012 Online only Revised for Version 6.17 (Release 2012a)
September 2012 Online only Revised for Version 6.18 (Release 2012b)
Contents
Filtering, Linear Systems and Transforms
Overview
1
Filter Implementation and Analysis ................. 1-2
Filtering Overview ................................ 1-2
Convolution and Filtering ........................... 1-2
Filters and Transfer Functions ...................... 1-3
Filtering with the filter Function ..................... 1-4
The filter Function ................................ 1-6
Other Functions for Filtering ....................... 1-8
Multirate Filter Bank Implementation ................ 1-8
Anti-Causal, Zero-Phase Filter Implementation ......... 1-9
Frequency Domain Filter Implementation ............. 1-10
Impulse Response ................................. 1-12
Frequency Response ............................... 1-14
DigitalDomain ................................... 1-14
Analog Domain ................................... 1-16
Magnitude and Phase .............................. 1-17
Delay ........................................... 1-18
Zero-Pole Analysis ................................. 1-21
Linear System Models .............................. 1-23
Available Models .................................. 1-23
Discrete-Time System Models ....................... 1-23
Continuous-Time System Models ..................... 1-31
Linear System Transformations ..................... 1-32
Discrete Fourier Transform ........................ 1-34
v
Filter Design and Implementation
2
Filter Requirements and Specification .............. 2-2
IIR Filter Design .................................. 2-4
IIR vs. FIR Filters ................................. 2-4
Classical IIR Filters ............................... 2-4
Other IIR Filters .................................. 2-5
IIR Filter Method Summary ......................... 2-5
Classical IIR Filter Design Using Analog Prototyping .... 2-6
Comparison of Classical IIR Filter Types .............. 2-9
FIR Filter Design .................................. 2-17
FIR vs. IIR Filters ................................. 2-17
FIR Filter Summary ............................... 2-18
Linear Phase Filters ............................... 2-18
Windowing Method ................................ 2-20
Multiband FIR Filter Design with Transition Bands ..... 2-24
Constrained Least Squares FIR Filter Design .......... 2-31
Arbitrary-Response Filter Design .................... 2-37
Special Topics in IIR Filter Design .................. 2-43
Classic IIR Filter Design ........................... 2-43
Analog Prototype Design ........................... 2-44
Frequency Transformation .......................... 2-44
Filter Discretization ............................... 2-46
Filtering Data With Signal Processing Toolbox
Software ........................................ 2-52
Lowpass FIR Filter — Window Method ................ 2-52
Bandpass Filters — Minimum-Order FIR and IIR
Systems ....................................... 2-56
Zero-Phase Filtering ............................... 2-65
Selected Bibliography .............................. 2-71
vi Contents
Designing a Filter in Fdesign — Process
Overview
3
Process Flow Diagram and Filter Design
Methodology .................................... 3-2
Exploring the Process Flow Diagram .................. 3-2
Selecting a Response ............................... 3-4
Selecting a Specification ............................ 3-4
Selecting an Algorithm ............................. 3-6
Customizing the Algorithm ......................... 3-8
Designing the Filter ............................... 3-8
Design Analysis ................................... 3-9
RealizeorApplytheFiltertoInputData .............. 3-10
Designing a Filter in the Filterbuilder GUI
4
Filterbuilder Design Process ....................... 4-2
Introduction to Filterbuilder ........................ 4-2
Design a Filter Using Filterbuilder ................... 4-2
Select a Response ................................. 4-3
Select a Specification .............................. 4-5
Select an Algorithm ................................ 4-5
Customize the Algorithm ........................... 4-6
Analyze the Design ................................ 4-8
RealizeorApplytheFiltertoInputData .............. 4-8
Designing a FIR Filter Using filterbuilder ........... 4-10
FIR Filter Design ................................. 4-10
FDATool: A Filter Design and Analysis GUI
5
Overview ......................................... 5-2
vii
FDATool ......................................... 5-2
Filter Design Methods ............................. 5-2
Using the Filter Design and Analysis Tool ............. 5-4
Analyzing Filter Responses ......................... 5-4
Filter Design and Analysis Tool Panels ................ 5-4
Getting Help ..................................... 5-5
Using FDATool .................................... 5-6
Choosing a Response Type .......................... 5-7
Choosing a Filter Design Method ..................... 5-8
Setting the Filter Design Specifications ............... 5-8
Computing the Filter Coefficients .................... 5-12
Analyzing the Filter ............................... 5-13
Editing the Filter Using the Pole/Zero Editor ........... 5-19
Converting the Filter Structure ...................... 5-23
Exporting a Filter Design ........................... 5-26
Generating a C Header File ......................... 5-32
Generating MATLAB Code .......................... 5-34
Managing Filters in the Current Session .............. 5-35
Saving and Opening Filter Design Sessions ............ 5-38
Importing a Filter Design .......................... 5-39
Import Filter Panel ................................ 5-39
Filter Structures .................................. 5-40
Statistical Signal Processing
6
Correlation and Covariance ........................ 6-2
Background Information ............................ 6-2
Using xcorr and xcov Functions ...................... 6-3
Bias and Normalization ............................ 6-3
Multiple Channels ................................. 6-4
Spectral Analysis .................................. 6-5
Background Information ............................ 6-5
Spectral Estimation Method ......................... 6-7
Nonparametric Methods ............................ 6-9
Parametric Methods ............................... 6-32
viii Contents
Selected Bibliography .............................. 6-46
Special Topics
7
Windows .......................................... 7-2
Why Use Windows? ................................ 7-2
Available Window Functions ........................ 7-2
Graphical User Interface Tools ...................... 7-3
Basic Shapes ..................................... 7-3
Generalized Cosine Windows ........................ 7-6
Kaiser Window ................................... 7-8
Chebyshev Window ................................ 7-12
Parametric Modeling .............................. 7-13
What is Parametric Modeling ........................ 7-13
Available Parametric Modeling Functions ............. 7-13
Time-Domain Based Modeling ....................... 7-14
Frequency-Domain Based Modeling .................. 7-18
Resampling ....................................... 7-21
Available Resampling Functions ..................... 7-21
resample Function ................................. 7-21
decimate and interp Functions ....................... 7-23
upfirdn Function .................................. 7-23
spline Function ................................... 7-23
Cepstrum Analysis ................................. 7-24
What Is a Cepstrum? .............................. 7-24
Inverse Complex Cepstrum ......................... 7-27
FFT-Based Time-Frequency Analysis ................ 7-28
Median Filtering .................................. 7-29
Communications Applications ...................... 7-30
Modulation ....................................... 7-30
Demodulation .................................... 7-31
ix
Voltage Controlled Oscillator ........................ 7-34
Deconvolution ..................................... 7-35
Specialized Transforms ............................ 7-36
Chirp z-Transform ................................. 7-36
Discrete Cosine Transform .......................... 7-37
Hilbert Transform ................................. 7-40
Walsh–Hadamard Transform ........................ 7-41
Selected Bibliography .............................. 7-47
SPTool: A Signal Processing GUI Suite
8
SPTool: An Interactive Signal Processing
Environment .................................... 8-2
SPTool Overview .................................. 8-2
SPTool Data Structures ............................ 8-3
Opening SPTool ................................... 8-4
Getting Context-Sensitive Help ..................... 8-6
Signal Browser .................................... 8-7
Overview of the Signal Browser ...................... 8-7
Opening the Signal Browser ......................... 8-7
FDATool .......................................... 8-10
Filter Visualization Tool ........................... 8-11
Connection between FVTool and SPTool ............... 8-11
Opening the Filter Visualization Tool ................. 8-11
Analysis Parameters ............................... 8-12
Spectrum Viewer .................................. 8-13
xContents
Spectrum Viewer Overview ......................... 8-13
Opening the Spectrum Viewer ....................... 8-13
Filtering and Analysis of Noise ..................... 8-16
Overview ........................................ 8-16
Importing a Signal into SPTool ...................... 8-16
Designing a Filter ................................. 8-18
Applying a Filter to a Signal ........................ 8-20
Analyzing a Signal ................................ 8-22
Spectral Analysis in the Spectrum Viewer ............. 8-24
Exporting Signals, Filters, and Spectra .............. 8-27
Opening the Export Dialog Box ...................... 8-27
Exporting a Filter to the MATLAB Workspace .......... 8-28
Accessing Filter Parameters ........................ 8-29
Accessing Filter Parameters in a Saved Filter .......... 8-29
Accessing Parameters in a Saved Spectrum ............ 8-30
Importing Filters and Spectra ...................... 8-32
Similarities to Other Procedures ..................... 8-32
Importing Filters .................................. 8-32
Importing Spectra ................................. 8-34
Loading Variables from the Disk .................... 8-36
Saving and Loading Sessions ....................... 8-37
SPTool Sessions ................................... 8-37
Filter Formats .................................... 8-37
Selecting Signals, Filters, and Spectra ............... 8-39
Editing Signals, Filters, or Spectra .................. 8-40
Making Signal Measurements with Markers ......... 8-41
Setting Preferences ................................ 8-43
Overview of Setting Preferences ..................... 8-43
Summary of Settable Preferences .................... 8-44
xi
Setting the Filter Design Tool ....................... 8-45
Using the Filter Designer ........................... 8-47
Filter Designer ................................... 8-47
Filter Types ...................................... 8-47
FIR Filter Methods ................................ 8-48
IIR Filter Methods ................................ 8-48
Pole/Zero Editor ................................... 8-48
Spectral Overlay Feature ........................... 8-48
Opening the Filter Designer ......................... 8-48
Accessing Filter Parameters in a Saved Filter .......... 8-50
Designing a Filter with the Pole/Zero Editor ........... 8-53
Positioning Poles and Zeros ......................... 8-54
Redesigning a Filter Using the Magnitude Plot ......... 8-56
Code Generation from MATLAB Support in
Signal Processing Toolbox
9
Supported Functions .............................. 9-2
Specifying Inputs in Code Generation from MATLAB
................................................. 9-7
Defining Input Size and Type ........................ 9-7
Inputs must be Constants ........................... 9-8
Code Generation Examples ......................... 9-11
Apply Window to Input Signal ....................... 9-11
Apply Lowpass Filter to Input Signal ................. 9-13
Cross Correlate or Autocorrelate Input Data ........... 9-14
freqz With No Output Arguments ................... 9-15
Zero Phase Filtering ............................... 9-16
xii Contents
Convolution and Correlation
10
Linear and Circular Convolution ................... 10-2
Confidence Intervals for Sample Autocorrelation ..... 10-5
Residual Analysis with Autocorrelation ............. 10-9
Autocorrelation of Moving Average Process .......... 10-19
Cross-correlation of Two Moving Average Processes .. 10-21
Cross-correlation of Delayed Signal in Noise ......... 10-23
Cross-correlation of Phase-Lagged Sine Wave ........ 10-26
Multirate Signal Processing
11
Downsampling — Signal Phases .................... 11-2
Downsampling — Aliasing .......................... 11-6
Filtering Before Downsampling ..................... 11-12
Upsampling — Imaging Artifacts .................... 11-14
Filtering After Upsampling — Interpolation ......... 11-16
Changing Signal Sampling Rate ..................... 11-18
xiii
Spectral Analysis
12
Power Spectral Density Estimates Using FFT ........ 12-2
Bias and Variability in the Periodogram ............. 12-10
Cross Spectrum and Magnitude-Squared Coherence .. 12-20
Amplitude Estimation and Zero Padding ............ 12-24
Significance Testing for Periodic Component ........ 12-27
Frequency Estimation by Subspace Methods ......... 12-29
Frequency-Domain Linear Regression ............... 12-32
Linear Prediction
13
Prediction Polynomial ............................. 13-2
Formant Estimation with LPC Coefficients .......... 13-5
AR Order Selection with Partial Autocorrelation
Sequence ....................................... 13-9
Transforms
14
Complex Cepstrum — Fundamental Frequency
Estimation ..................................... 14-2
xiv Contents
Analytic Signal for Cosine .......................... 14-6
Envelope Extraction Using The Analytic Signal ...... 14-9
Signal Generation
15
Display Time-Domain Data in Signal Browser ........ 15-2
Import and Display Signals ......................... 15-3
Configure the Signal Browser Properties .............. 15-6
Modify the Signal Browser Display ................... 15-9
Inspect Your Data (Scaling the Axes and Zooming) ...... 15-11
Signal Measurement
16
RMS Value of Periodic Waveforms .................. 16-2
Slew Rate of Triangular Waveform .................. 16-5
Duty Cycle of Rectangular Pulse Waveform .......... 16-9
Estimate State for Digital Clock ..................... 16-12
CalculateSettlingTimewithSignalBrowser ......... 16-16
Find Peak Amplitudes in Signal Browser ............ 16-20
xv
Technical Conventions
A
Index
xvi Contents
1
Filtering, Linear Systems
and Transforms Overview
“Filter Implementation and Analysis” on page 1-2
“The filter Function” on page 1-6
“Other Functions for Filtering” on page 1-8
“Impulse Response” on page 1-12
“Frequency Response” on page 1-14
“Zero-Pole Analysis” on page 1-21
“Linear System Models” on page 1-23
“Discrete Fourier Transform” on page 1-34
1Filtering, Linear Systems and Transforms Overview
Filter Implementation and Analysis
In this section...
“Filtering Overview” on page 1-2
“Convolution and Filtering” on page 1-2
“Filters and Transfer Functions” on page 1-3
“Filtering with the filter Function” on page 1-4
Filtering Overview
This section describes how to filter discrete signals using the MATLAB®
filter function and other Signal Processing Toolbox™ functions. It also
discusses how to use the toolbox functions to analyze filter characteristics,
including impulse response, magnitude and phase response, group delay,
and zero-pole locations.
Convolution and Filtering
The mathematical foundation of filtering is convolution. The MATLAB conv
function performs standard one-dimensional convolution, convolving one
vector with another:
conv([1 1 1],[1 1 1])
ans =
12321
Note Convolve rectangular matrices for two-dimensional signal processing
using the conv2 function.
A digital filter’s output y(k)isrelatedtoitsinputx(k) by convolution with its
impulse response h(k).
yk hlxk l
l
() ()( )=−
=−
1-2
Filter Implementation and Analysis
If a digital filter’s impulse response h(k) is finite in length, and the input x(k)
is also of finite length, you can implement the filter using conv.Storex(k)ina
vector x,h(k)inavectorh, and convolve the two:
x = randn(5,1); % A random vector of length 5
h = [1 1 1 1]/4; % Length 4 averaging filter
y = conv(h,x);
The length of the output is the sum of the finite-length input vectors minus 1.
Filters and Transfer Functions
In general, the z-transform Y(z) of a discrete-time filter’s output y(n)isrelated
to the z-transform X(z)oftheinputby
Yz HzXz bbz bnz
aaz a
n
() () () () () ... ( )
() () ...
==
++++
+++
−−
12 1
12
1
1(() ()
mz
Xz
m
+
1
where H(z)isthefilterstransfer function. Here, the constants b(i)anda(i)are
the filter coefficients and the order of the filter is the maximum of nand m.
Note The filter coefficients start with subscript 1, rather than 0. This reflects
the standard indexing scheme used for MATLAB vectors.
MATLAB filter functions store the coefficients in two vectors, one for the
numerator and one for the denominator. By convention, it uses row vectors
forfiltercoefficients.
Filter Coefficients and Filter Names
Many standard names for filters reflect the number of aand bcoefficients
present:
When n=0(thatis,bis a scalar), the filter is an Infinite Impulse Response
(IIR), all-pole, recursive, or autoregressive (AR) filter.
When m=0(thatis,ais a scalar), the filter is a Finite Impulse Response
(FIR), all-zero, nonrecursive, or moving-average (MA) filter.
1-3
1Filtering, Linear Systems and Transforms Overview
If both nand mare greater than zero, the filter is an IIR, pole-zero,
recursive, or autoregressive moving-average (ARMA) filter.
The acronyms AR, MA, and ARMA are usually applied to filters associated
with filtered stochastic processes.
Filtering with the filter Function
It is simple to work back to a difference equation from the z-transform relation
shown earlier. Assume that a(1) = 1. Move the denominator to the left-hand
side and take the inverse z-transform.
yk a yk am yk m b xk b xk bn() ()( ) ( )( ) ()() ()( ) (+ − +…+ + = + − +…+ +21 1 1 21 1))( )xk n
In terms of current and past inputs, and past outputs, y(k)is
yk a ykb xk b xk bn xk n am() ()( )() () () ( ) ( ) ( ) (=++++− …+−−121 1 121 ))( )yk m
This is the standard time-domain representation of a digital filter, computed
starting with y(1) and assuming a causal system with zero initial conditions.
This representation’s progression is
ybx
ybxbxay
ybx
() () ()
() ()() ()() ()()
() ()()
111
2122121
313
=
=+
=+bbx bx ay ay()() ()() ()() ()()22 31 22 31+− −
=
A filter in this form is easy to implement with the filter function. For
example, a simple single-pole filter (lowpass) is
B = 1; % Numerator
A = [1 -0.9]; % Denominator
where the vectors Band Arepresent the coefficients of a filter in transfer
function form. Note that the Acoefficient vectors are written as if the output
and input terms are separated in the difference equation. For the example,
the previous coefficient vectors represent a linear constant-coefficient
difference equation of
1-4
Filter Implementation and Analysis
yn yn xn() . ( ) ()−−=09 1
Changing the sign of the A(2) coefficient, results in the difference equation
yn yn xn() . ( ) ()+−=09 1
The previous coefficients are represented as:
B = 1; %Numerator
A = [1 0.9]; %Denominator
and results in a highpass filter.
Toapplythisfiltertoyourdata,use
y = filter(B,A,x);
filter gives you as many output samples as there are input samples, that
is, the length of yis the same as the length of x.Ifthefirstelementofa
is not 1, filter divides the coefficients by a(1) before implementing the
difference equation.
1-5
1Filtering, Linear Systems and Transforms Overview
The filter Function
filter is implemented as the transposed direct-form II structure, where n-1
is the filter order. This is a canonical form that has the minimum number of
delay elements.
At sample m,filter computes the difference equations
ym b xm z m
zm b xm zm a ym
zn
() ()() ( )
() ()() ( ) ()()
=+
=+
=
11
212
1
12

221
1
111() ( )() ( ) ( )()
() ()() ()
mbn xmz m an ym
zmbnxman
n
n
=− + −−
=−
yym()
In its most basic form, filter initializes the delay outputs zi(1), i= 1, ..., n-1
to 0. This is equivalent to assuming both past inputs and outputs are zero.
Set the initial delay outputs using a fourth input parameter to filter,or
access the final delay outputs using a second output parameter:
[y,zf] = filter(b,a,x,zi)
Access to initial and final conditions is useful for filtering data in sections,
especially if memory limitations are a consideration. Suppose you have
collected data in two segments of 5000 points each:
x1 = randn(5000,1); % Generate two random data sequences.
x2 = randn(5000,1);
Perhaps the first sequence, x1, corresponds to the first 10 minutes of data
and the second, x2, to an additional 10 minutes. The whole sequence is
1-6
The filter Function
x=[x1;x2]. If there is not sufficient memory to hold the combined sequence,
filter the subsequences x1 and x2 one at a time. To ensure continuity of
the filtered sequences, use the final conditions from x1 as initial conditions
to filter x2:
[y1,zf] = filter(b,a,x1);
y2 = filter(b,a,x2,zf);
The filtic function generates initial conditions for filter.filtic computes
the delay vector to make the behavior of the filter reflect past inputs and
outputs that you specify. To obtain the same output delay values zf as above
using filtic,use
zf = filtic(b,a,flipud(y1),flipud(x1));
This can be useful when filtering short data sequences, as appropriate initial
conditions help reduce transient startup effects.
1-7
1Filtering, Linear Systems and Transforms Overview
Other Functions for Filtering
In this section...
“Multirate Filter Bank Implementation” on page 1-8
“Anti-Causal, Zero-Phase Filter Implementation” on page 1-9
“Frequency Domain Filter Implementation” on page 1-10
Multirate Filter Bank Implementation
The upfirdn function alters the sampling rate of a signal by an integer ratio
P/Q. It computes the result of a cascade of three systems that performs the
following tasks:
Upsampling (zero insertion) by integer factor p
FilteringbyFIRfilterh
Downsampling by integer factor q
For example, to change the sample rate of a signal from 44.1 kHz to 48 kHz,
we first find the smallest integer conversion ratio p/q.Set
d = gcd(48000,44100);
p = 48000/d;
q = 44100/d;
In this example, p=160 and q=147. Sample rate conversion is then
accomplished by typing
y = upfirdn(x,h,p,q)
This cascade of operations is implemented in an efficient manner using
polyphase filtering techniques, and it is a central concept of multirate
filtering. Note that the quality of the resampling result relies on the quality of
the FIR filter h.
1-8
Other Functions for Filtering
Filter banks may be implemented using upfirdn by allowing the filter h
to be a matrix, with one FIR filter per column. A signal vector is passed
independently through each FIR filter, resulting in a matrix of output signals.
Other functions that perform multirate filtering (with fixed filter) include
resample,interp,anddecimate.
Anti-Causal, Zero-Phase Filter Implementation
In the case of FIR filters, it is possible to design linear phase filters that,
when applied to data (using filter or conv), simply delay the output by a
fixed number of samples. For IIR filters, however, the phase distortion is
usually highly nonlinear. The filtfilt function uses the information in the
signal at points before and after the current point, in essence “looking into the
future,” to eliminate phase distortion.
To see how filtfilt does this, recall that if the z-transform of a real
sequence x(n)isX(z), the z-transform of the time reversed sequence x(n)is
X(1/z). Consider the processing scheme.
Image of Anti Causal Zero Phase Filter
When |z|=1,thatisz=ejω, the output reduces to X(ejω)|H(ejω)|2.Given
all the samples of the sequence x(n), a doubly filtered version of xthat has
zero-phase distortion is possible.
For example, a 1-second duration signal sampled at 100 Hz, composed of two
sinusoidal components at 3 Hz and 40 Hz, is
fs = 100;
t = 0:1/fs:1;
x = sin(2*pi*t*3)+.25*sin(2*pi*t*40);
Now create a 10-point averaging FIR filter, and filter xusing both filter
and filtfilt for comparison:
1-9
1Filtering, Linear Systems and Transforms Overview
b = ones(1,10)/10; % 10 point averaging filter
y = filtfilt(b,1,x); % Noncausal filtering
yy = filter(b,1,x); % Normal filtering
plot(t,x,t,y,'--',t,yy,':')
Both filtered versions eliminate the 40 Hz sinusoid evident in the original,
solid line. The plot also shows how filter and filtfilt differ; the dashed
(filtfilt) line is in phase with the original 3 Hz sinusoid, while the dotted
(filter) line is delayed by about five samples. Also, the amplitude of the
dashed line is smaller due to the magnitude squared effects of filtfilt.
filtfilt reduces filter startup transients by carefully choosing initial
conditions, and by prepending onto the input sequence a short, reflected piece
of the input sequence. For best results, make sure the sequence you are
filtering has length at least three times the filter order and tapers to zero on
both edges.
Frequency Domain Filter Implementation
Dualitybetweenthetimedomainandthefrequencydomainmakesitpossible
to perform any operation in either domain. Usually one domain or the other is
more convenient for a particular operation, but you can always accomplish
a given operation in either domain.
1-10
Other Functions for Filtering
To implement general IIR filtering in the frequency domain, multiply the
discrete Fourier transform (DFT) of the input sequence with the quotient of
the DFT of the filter:
n = length(x);
y = ifft(fft(x).*fft(b,n)./fft(a,n));
This computes results that are identical to filter, but with different startup
transients (edge effects). For long sequences, this computation is very
inefficient because of the large zero-padded FFT operations on the filter
coefficients, and because the FFT algorithm becomes less efficient as the
number of points nincreases.
For FIR filters, however, it is possible to break longer sequences into shorter,
computationally efficient FFT lengths. The function
y = fftfilt(b,x)
uses the overlap add method to filter a long sequence with multiple
medium-length FFTs. Its output is equivalent to filter(b,1,x).
1-11
1Filtering, Linear Systems and Transforms Overview
Impulse Response
The impulse response of a digital filter is the output arising from the unit
impulse sequence defined as
()nn
n
==
10
00
You can generate an impulse sequence a number of ways; one straightforward
way is
imp = [1; zeros(49,1)];
The impulse response of the simple filter b=1and a=[1 -0.9] is
h = filter(b,a,imp);
Asimplew
ay to display the impulse response is with the Filter Visualization
Tool (fvtool):
fvtool(b,a)
Then click the Impulse Response button on the toolbar or select
Analysis > Impulse Response. This plot shows the exponential decay
h(n) =0.9n ofthesinglepolesystem:
1-12
Impulse Response
1-13
1Filtering, Linear Systems and Transforms Overview
Frequency Response
In this section...
“Digital Domain” on page 1-14
“Analog Domain” on page 1-16
“Magnitude and Phase” on page 1-17
“Delay” on page 1-18
Digital Domain
freqz uses an FFT-based algorithm to calculate the z-transform frequency
response of a digital filter. Specifically, the statement
[h,w] = freqz(b,a,p)
returns the p-point complex frequency response, H(ejω), of the digital filter.
He bbe bne
aae am
jw jw jwn
jw
() ( ) ( ) ... ( )
() () ... (
=++++
++++
−−
12 1
12 11)ejwm
In its simplest form, freqz accepts the filter coefficient vectors band a,andan
integer pspecifying the number of points at which to calculate the frequency
response. freqz returns the complex frequency response in vector h,andthe
actual frequency points in vector win rad/s.
freqz can accept other parameters, such as a sampling frequency or a
vector of arbitrary frequency points. The example below finds the 256-point
frequency response for a 12th-order Chebyshev Type I filter. The call to freqz
specifies a sampling frequency fs of 1000 Hz:
[b,a] = cheby1(12,0.5,200/500);
[h,f] = freqz(b,a,256,1000);
Because the parameter list includes a sampling frequency, freqz returns a
vector fthat contains the 256 frequency points between 0 and fs/2 used in
the frequency response calculation.
1-14
Frequency Response
Note This toolbox uses the convention that unit frequency is the Nyquist
frequency, defined as half the sampling frequency. The cutoff frequency
parameter for all basic filter design functions is normalized by the Nyquist
frequency. For a system with a 1000 Hz sampling frequency, for example,
300 Hz is 300/500 = 0.6. To convert normalized frequency to angular frequency
around the unit circle, multiply by π. To convert normalized frequency back
to hertz, multiply by half the sample frequency.
If you call freqz with no output arguments, it plots both magnitude
versus frequency and phase versus frequency. For example, a ninth-order
Butterworth lowpass filter with a cutoff frequency of 400 Hz, based on a 2000
Hz sampling frequency, is
[b,a] = butter(9,400/1000);
To calculate the 256-point complex frequency response for this filter, and plot
the magnitude and phase with freqz,use
freqz(b,a,256,2000)
or to display the magnitude and phase responses in fvtool, which provides
additional analysis tools, use
fvtool(b,a)
and click the Magnitude and Phase Response button on the toolbar or
select Analysis > Magnitude and Phase Response.
1-15
1Filtering, Linear Systems and Transforms Overview
freqz can also accept a vector of arbitrary frequency points for use in the
frequency response calculation. For example,
w = linspace(0,pi);
h = freqz(b,a,w);
calculates the complex frequency response at the frequency points in wfor the
filter defined by vectors band a. The frequency points can range from 0 to
2π. To specify a frequency vector that ranges from zero to your sampling
frequency, include both the frequency vector and the sampling frequency
value in the parameter list.
Analog Domain
freqs evaluates frequency response for an analog filter defined by two input
coefficient vectors, band a. Its operation is similar to that of freqz;youcan
specify a number of frequency points to use, supply a vector of arbitrary
frequency points, and plot the magnitude and phase response of the filter.
1-16
Frequency Response
Magnitude and Phase
MATLAB functions are available to extract magnitude and phase from a
frequency response vector h. The function abs returns the magnitude of the
response; anglereturns the phase angle in radians. To extract the magnitude
and phase of a Butterworth filter:
[b,a] = butter(9,400/1000);
fvtool(b,a)
and click the Magnitude and Phase Response button on the toolbar or
select Analysis > Magnitude and Phase Response to display the plot.
The unwrap function is also useful in frequency analysis. unwrap unwraps
the phase to make it continuous across 360º phase discontinuities by adding
multiples of ±360°, as needed. To see how unwrap is useful, design a
25th-order lowpass FIR filter:
h = fir1(25,0.4);
Obtain the filter’s frequency response with freqz, and plot the phase in
degrees:
1-17
1Filtering, Linear Systems and Transforms Overview
[H,f] = freqz(h,1,512,2);
plot(f,angle(H)*180/pi); grid
It is difficult to distinguish the 360° jumps (an artifact of the arctangent
function inside angle) from the 180° jumps that signify zeros in the frequency
response.
unwrap eliminates the 360° jumps:
plot(f,unwrap(angle(H))*180/pi);
or you can use phasez to see the unwrapped phase.
Delay
The group delay ofafilterisameasureoftheaveragetimedelayofthefilter
as a function of frequency. It is defined as the negative first derivative of a
filter’s phase response. If the complex frequency response of a filter is H(ejω),
then the group delay is
 
g
d
d
() ()
=−
1-18
Frequency Response
where θ(ω) is the phase, or argument of H(ejω). Compute group delay with
[gd,w] = grpdelay(b,a,n)
which returns the n-point group delay, ,τg(ω) of the digital filter specified by
band a, evaluated at the frequencies in vector w.
The phase delay of a filter is the negative of phase divided by frequency:
 
p() ()
=−
To plot both the group and phase delays of a system on the same FVTool
graph, type
[b,a] = butter(10,200/1000);
hFVT = fvtool(b,a,'Analysis','grpdelay');
set(hFVT,'NumberofPoints',128,'OverlayedAnalysis','phasedelay');
legend(hFVT)
1-19
1Filtering, Linear Systems and Transforms Overview
1-20
Zero-Pole Analysis
Zero-Pole Analysis
The zplane function plots poles and zeros of alinearsystem. Forexample,
a simple filter with a zero at -1/2 and a complex pole pair at 0.9e–j2π(0.3) and
0.9ej2π(0.3) is
zer = -0.5;
pol = 0.9*exp(j*2*pi*[-0.3 0.3]');
To view the pole-zero plot for this filter you can use
zplane(zer,pol)
or, for access to additional tools, use fvtool. First convert the poles and zeros
to transfer function form, then call fvtool,
[b,a] = zp2tf(zer,pol,1);
fvtool(b,a)
and click the Pole/Zero Plot toolbar button on the toolbar or select
Analysis > Pole/Zero Plot to see the plot.
1-21
1Filtering, Linear Systems and Transforms Overview
For a system in zero-pole form, supply column vector arguments zand pto
zplane:
zplane(z,p)
For a system in transfer function form, supply row vectors band aas
arguments to zplane:
zplane(b,a)
In this case zplane finds the roots of band ausing the roots function and
plots the resulting zeros and poles.
See “Linear System Models” on page 1-23 for details on zero-pole and transfer
function representation of systems.
1-22
Linear System Models
Linear System Models
In this section...
“Available Models” on page 1-23
“Discrete-Time System Models” on page 1-23
“Continuous-Time System Models” on page 1-31
“Linear System Transformations” on page 1-32
Available Models
Several Signal Processing Toolbox models are provided for representing linear
time-invariant systems. This flexibility lets you choose the representational
scheme that best suits your application and, within the bounds of numeric
stability, convert freely to and from most other models. This section provides
a brief overview of supported linear system models and describes how to work
with these models in the MATLAB technical computing environment.
Discrete-Time System Models
The discrete-time system models are representational schemes for digital
filters. The MATLAB technical computing environment supports several
discrete-time system models, which are described in the following sections:
“Transfer Function” on page 1-23
“Zero-Pole-Gain” on page 1-24
“State-Space” on page 1-25
“Partial Fraction Expansion (Residue Form)” on page 1-26
“Second-Order Sections (SOS)” on page 1-27
“Lattice Structure” on page 1-28
“Convolution Matrix” on page 1-30
Transfer Function
The transfer function is a basic z-domain representation of a digital filter,
expressing the filter as a ratio of two polynomials. It is the principal
1-23
1Filtering, Linear Systems and Transforms Overview
discrete-time model for this toolbox. The transfer function model description
for the z-transform of a digital filter’s difference equation is
Yz bbz bnz
aaz amz
Xz
n
m
() () () ( )
() () ( ) ()=++++
++++
−−
−−
12 1
12 1
1
1
Here, the constants b(i)anda(i) are the filter coefficients, and the order of the
filter is the maximum of nand m. In the MATLAB environment, you store
these coefficients in two vectors (row vectors by convention), one row vector
for the numerator and one for the denominator. See “Filters and Transfer
Functions” on page 1-3 for more details on the transfer function form.
Zero-Pole-Gain
The factored or zero-pole-gain form of a transfer function is
Hz qz
pz kzq zq zqn
zp zp
() ()
()
( ( ))( ( ))...( ( ))
( ( ))( ( ))
==
−− −
−−
12
12
....( ( ))zpn
By convention, polynomial coefficients are stored in row vectors and
polynomial roots in column vectors. In zero-pole-gain form, therefore, the
zero and pole locations for the numerator and denominator of a transfer
function reside in column vectors. The factored transfer function gain kis a
MATLAB scalar.
The poly and roots functions convert between polynomial and zero-pole-gain
representations. For example, a simple IIR filter is
b=[234];
a=[1331];
The zeros and poles of this filter are
q = roots(b)
p = roots(a)
% Gain factor
k = b(1)/a(1)
Returning to the original polynomials,
1-24
Linear System Models
bb = k*poly(q)
aa = poly(p)
Note that band ain this case represent the transfer function:
Hz zz
zzz
zz
zzz
()=++
+++
=++
+++
−−
−−
23 4
13 3
234
331
12
123
2
32
For b=[2 3 4],theroots function misses the zero for zequal to 0. In fact, it
misses poles and zeros for zequal to 0 whenever the input transfer function
has more poles than zeros, or vice versa. This is acceptable in most cases. To
circumvent the problem, however, simplyappendzerostomakethevectors
the same length before using the roots function; for example, b=[b 0].
State-Space
It is always possible to represent a digital filter, or a system of difference
equations, as a set of first-order difference equations. In matrix or state-space
form, you can write the equations as
xn Axn Bun
yn Cxn Dun
( ) () ()
() () ()
+= +
=+
1
where uis the input, xis the state vector, and yis the output. For
single-channel systems, Ais an m-by-mmatrix where mis the order of the filter,
Bis a column vector, Cis a row vector, and Dis a scalar. State-space notation
is especially convenient for multichannel systems where input uand output y
become vectors, and B,C,andDbecome matrices.
State-space representation extends easily to the MATLAB environment.A,B,
C,andDare rectangular arrays; MATLAB functions treat them as individual
variables.
Taking the z-transform of the state-space equations and combining them
shows the equivalence of state-space and transfer function forms:
Yz HzUz Hz CzI A B D() () () () ( )==+
, where 1
1-25
1Filtering, Linear Systems and Transforms Overview
Don’t be concerned if you are not familiar with the state-space representation
of linear systems. Some of the filter design algorithms use state-space
form internally but do not require any knowledge of state-space concepts
to use them successfully. If your applications use state-space based signal
processing extensively, however, see the Control System Toolbox™ product
for a comprehensive library of state-space tools.
Partial Fraction Expansion (Residue Form)
Each transfer function also has a corresponding partial fraction expansion or
residue form representation, given by
bz
az
r
pz
rn
pnz
kkz
()
()
()
() ... ()
() ( ) ( ) ...=++
++ ++
−−
1
11 1 12
11
1kkm n z mn
()
()
−+ −−
1
provided H(z) has no repeated poles. Here, nis the degree of the denominator
polynomial of the rational transfer function b(z)/a(z). If ris a pole of
multiplicity sr,thenH(z)hastermsoftheform:
rj
pjz
rj
pjz
rj s
pjz
r
sr
()
()
()
(())
... ()
(())1
1
1
1
1
112 1
++
++−
−− −
The Signal Processing Toolbox residuez function in converts transfer
functions to and from the partial fraction expansion form. The “z”ontheend
of residuez stands for z-domain, or discrete domain. residuez returns the
poles in a column vector p, the residues corresponding to the poles in a column
vector r, and any improper part of the original transfer function in a row
vector k.residuez determines that two poles are the same if the magnitude of
their difference is smaller than 0.1 percent of either of the poles’ magnitudes.
Partial fraction expansion arises in signal processing as one method of finding
the inverse z-transform of a transfer function. For example, the partial
fraction expansion of
Hz z
zz
()=−+
++
−−
48
16 8
1
12
is
1-26
Linear System Models
b = [-4 8];
a=[168];
[r,p,k] = residuez(b,a)
which corresponds to
Hz
zz
()=
+++
−−
12
14
8
12
11
To find the inverse z-transform of H(z), find the sum of the inverse
z-transforms of the two addends of H(z), giving the causal impulse response:
hn n
nn
() ( ) ( ) ,,,=− + − = 12 4 8 2 0 1 2
To verify this in the MATLAB environment, type
imp=[10000];
resptf = filter(b,a,imp)
respres = filter(r(1),[1 -p(1)],imp)+...
filter(r(2),[1 -p(2)],imp)
Second-Order Sections (SOS)
Any transfer function H(z) has a second-order sections representation
Hz H z bbzbz
aazaz
k
k
Lkk k
kk k
k
L
() ()==
++
++
=
−−
−−
=
∏∏
1
01
122
01
122
1
where Lis the number of second-order sections that describe the system.
The MATLAB environment represents the second-order section form of a
discrete-time system as an L-by-6 array sos. Each row of sos contains a
single second-order section, where the row elements are the three numerator
and three denominator coefficients that describe the second-order section.
1-27
1Filtering, Linear Systems and Transforms Overview
sos
bbbaaa
bbbaaa
bb
LL
=
01 11 21 01 11 21
02 12 22 02 12 22
01
......
......
bbaaa
LLLL2012
There are many ways to represent a filter in second-order section form.
Through careful pairing of the pole and zero pairs, ordering of the sections
in the cascade, and multiplicative scaling of the sections, it is possible to
reduce quantization noise gain and avoid overflow in some fixed-point filter
implementations. The functions zp2sos and ss2sos, described in “Linear
System Transformations” on page 1-32, perform pole-zero pairing, section
scaling, and section ordering.
Note AllSignalProcessingToolboxsecond-order section transformations
apply only to digital filters.
Lattice Structure
For a discrete Nth order all-pole or all-zero filter described by the polynomial
coefficients a(n),n=1,2,...,N+1, there are Ncorresponding lattice structure
coefficients k(n),n=1,2,...,N. The parameters k(n) are also called the
reflection coefficients of the filter. Given these reflection coefficients, you can
implement a discrete filter as shown below.
1-28
Linear System Models
FIRandIIRLatticeFilterstructurediagrams
For a general pole-zero IIR filter described by polynomial coefficients a
and b, there are both lattice coefficients k(n) for the denominator aand
ladder coefficients v(n) for the numerator b. The lattice/ladder filter may be
implemented as
Diagram of lattice/ladder filter
The toolbox function tf2latc accepts an FIR or IIR filter in polynomial form
and returns the corresponding reflection coefficients. An example FIR filter
in polynomial form is
b = [1.0000 0.6149 0.9899 0.0000 0.0031 -0.0082];
1-29
1Filtering, Linear Systems and Transforms Overview
This filter’s lattice (reflection coefficient) representation is
k = tf2latc(b)
For IIR filters, the magnitude of the reflection coefficients provides an easy
stability check. If all the reflection coefficients corresponding to a polynomial
have magnitude less than 1, all of that polynomial’s roots are inside the unit
circle. For example, consider an IIR filter with numerator polynomial bfrom
above and denominator polynomial:
a = [1 1/2 1/3];
The filter’s lattice representation is
[k,v] = tf2latc(b,a);
Because abs(k) <1for all reflection coefficients in k, the filter is stable.
The function latc2tf calculates the polynomial coefficients for a filter from
its lattice (reflection) coefficients. Given the reflection coefficient vector
k(above), the corresponding polynomial form is
b = latc2tf(k);
The lattice or lattice/ladder coefficients can be used to implement the filter
using the function latcfilt.
Convolution Matrix
In signal processing, convolving two vectors or matrices is equivalent to
filtering one of the input operands by the other. This relationship permits the
representation of a digital filter as a convolution matrix.
Given any vector, the toolbox function convmtx generates a matrix whose
inner product with another vector is equivalent to the convolution of the two
vectors. The generated matrix represents a digital filter that you can apply to
any vector of appropriate length; the inner dimension of the operands must
agreetocomputetheinnerproduct.
The convolution matrix for a vector b, representing the numerator coefficients
for a digital filter, is
1-30
Linear System Models
b = [1 2 3]; x = randn(3,1);
C = convmtx(b',3);
Two equivalent ways to convolve bwith xare as follows.
y1 = C*x;
y2 = conv(b,x);
Continuous-Time System Models
The continuous-time system models are representational schemes for analog
filters. Many of the discrete-time system models described earlier are also
appropriate for the representation of continuous-time systems:
State-space form
Partial fraction expansion
Transfer function
Zero-pole-gain form
It is possible to represent any system of linear time-invariant differential
equations as a set of first-order differential equations. In matrix or state-space
form, you can express the equations as
xAxBu
yCxDu
=+
=+
where uis a vector of nu inputs, xis an nx-element state vector, and yis a
vector of ny outputs. In the MATLAB environment, A,B,C,andDare stored in
separate rectangular arrays.
An equivalent representation of the state-space system is the Laplace
transform transfer function description
Ys HsUs() () ()=
where
1-31
1Filtering, Linear Systems and Transforms Overview
Hs CsI A B D() ( )=− +
1
For single-input, single-output systems, this form is given by
Hs bs
as
bs bs bn
as as a
nn
mm
() ()
()
() () ( )
() () (
== ++++
+++
12 1
12
1
1mm +1)
Given the coefficients of a Laplace transform transfer function, residue
determines the partial fraction expansion of the system. See the description
of residue for details.
The factored zero-pole-gain form is
Hs zs
ps ksz sz szn
sp sp
() ()
()
( ( ))( ( )) ( ( ))
( ( ))( ( )) (
==−−
−−
12
12
sspm())
Asinthediscrete-timecase,theMATLABenvironmentstorespolynomial
coefficients in row vectors in descending powers of s.Itstorespolynomial
roots, or zeros and poles, in column vectors.
Linear System Transformations
A number of Signal Processing Toolbox functions are provided to convert
between the various linear system models.. You can use the following chart to
find an appropriate transfer function: find the row of the model to convert
from on the left side of the chart and the column of the model to convert to on
the top of the chart and read the function name(s) at the intersection of the
row and column. Note that some cells of this table are empty.
Transfer
Function
State-
Space
Zero-
Pole-
Gain
Partial
Fraction
Lattice
Filter
Second-
Order
Sections
Convolution
Matrix
Transfer
Function
tf2ss tf2zp
roots
residuez tf2latc none convmtx
State-Space ss2tf ss2zp none none ss2sos none
1-32
Linear System Models
Transfer
Function
State-
Space
Zero-
Pole-
Gain
Partial
Fraction
Lattice
Filter
Second-
Order
Sections
Convolution
Matrix
Zero-Pole-
Gain
zp2tf
poly
zp2ss none none zp2sos none
Partial
Fraction
residuez none none none none none
Lattice
Filter
latc2tf none none none none none
SOS sos2tf sos2ss sos2zp none none none
Note Converting from one filter structure or model to another may produce
a result with different characteristics than the original. This is due to the
computer’s finite-precision arithmetic and the variations in the conversion’s
round-off computations.
Many of the toolbox filter design functions use these functions internally.
For example, the zp2ss function converts the poles and zeros of an analog
prototype into the state-space form required for creation of a Butterworth,
Chebyshev, or elliptic filter. Once in state-space form, the filter design
function performs any required frequency transformation, that is, it
transforms the initial lowpass design into a bandpass, highpass, or bandstop
filter, or a lowpass filter with the desired cutoff frequency.
Note AllSignalProcessingToolboxsecond-order section transformations
apply only to digital filters.
1-33
1Filtering, Linear Systems and Transforms Overview
Discrete Fourier Transform
The discrete Fourier transform, or DFT, is the primary tool of digital signal
processing. The foundation of Signal Processing Toolbox product is the fast
Fourier transform (FFT), a method for computing the DFT with reduced
execution time. Many of the toolbox functions (including z-domain frequency
response, spectrum and cepstrum analysis, and some filter design and
implementation functions) incorporate the FFT.
The MATLAB environment provides the functions fft and ifft to compute
the discrete Fourier transform and its inverse, respectively. For the input
sequence xand its transformed version X(the discrete-time Fourier transform
at equally spaced frequencies around the unit circle), the two functions
implement the relationships
Xk xn W
n
Nkn
N
() (),+= +
=
11
0
1
and
xn NXk W
k
N
Nkn
() .() 
111
0
1
In these equations, the series subscripts begin with 1 instead of 0 because of
the MATLAB vector indexing scheme, and
We
NjN
=2
/.
Note TheMATLABconventionistouseanegativejfor the fft function.
This is an engineering convention; physics and pure mathematics typically
use a positive j.
fft, with a single input argument x, computes the DFT of the input vector
or matrix. If xis a vector, fft computes the DFT of the vector; if xis a
rectangular array, fft computes the DFT of each array column.
1-34
Discrete Fourier Transform
For example, create a time vector and signal:
t = (0:1/100:10-1/100); % Time vector
x = sin(2*pi*15*t) + sin(2*pi*40*t); % Signal
The DFT of the signal, and the magnitude and phase of the transformed
sequence, are then
y = fft(x); % Compute DFT of x
m = abs(y); p = unwrap(angle(y)); % Magnitude and phase
To plot the magnitude and phase, type the following commands:
f = (0:length(y)-1)*99/length(y); % Frequency vector
plot(f,m); title('Magnitude');
set(gca,'XTick',[15 40 60 85]);
figure; plot(f,p*180/pi); title('Phase');
set(gca,'XTick',[15 40 60 85]);
Aseco
nd argument to fft specifies a number of points nfor the transform,
representing DFT length:
y = fft(x,n);
In this case, fft pads the input sequence with zeros if it is shorter than n,or
truncates the sequence if it is longer than n.Ifnisnotspecified,itdefaults
to the length of the input sequence. Execution time for fft depends on the
length, n, of the DFT it performs; see the fft for details about the algorithm.
1-35
1Filtering, Linear Systems and Transforms Overview
Note The resulting FFT amplitude is A*n/2, where A is the original amplitude
and n is the number of FFT points. This is true only if the number of FFT
points is greater than or equal to the number of data samples. If the number
of FFT points is less, the FFT amplitude is lower than the original amplitude
by the above amount.
The inverse discrete Fourier transform function ifft also accepts an input
sequence and, optionally, the number of desired points for the transform. Try
the example below; the original sequence xand the reconstructed sequence
are identical (within rounding error).
t = (0:1/255:1);
x = sin(2*pi*120*t);
y = real(ifft(fft(x)));
This toolbox also includes functions for the two-dimensional FFT and its
inverse, fft2 and ifft2. These functions are useful for two-dimensional
signal or image processing. The goertzel function, which is another
algorithm to compute the DFT, also is included in the toolbox. This function is
efficient for computing the DFTofaportionofalongsignal.
It is sometimes convenient to rearrange the output of the fft or fft2 function
so the zero frequency component is at the center of the sequence. The
MATLAB function fftshift moves the zero frequency component to the
center of a vector or matrix.
1-36
2
Filter Design and
Implementation
“Filter Requirements and Specification” on page 2-2
“IIR Filter Design” on page 2-4
“FIR Filter Design” on page 2-17
“Special Topics in IIR Filter Design” on page 2-43
“Filtering Data With Signal Processing Toolbox Software” on page 2-52
“Selected Bibliography” on page 2-71
2Filter Design and Implementation
Filter Requirements and Specification
Filter design is the process of creating the filter coefficients to meet specific
filtering requirements. Filter implementation involves choosing and applying
a particular filter structure to those coefficients. Only after both design and
implementation have been performed can data be filtered. The following
chapter describes filter design and implementation in Signal Processing
Toolbox software.
The goal of filter design is to perform frequency dependent alteration of a data
sequence. A possible requirement might be to remove noise above 200 Hz from
a data sequence sampled at 1000 Hz. A more rigorous specification might call
for a specific amount of passband ripple, stopband attenuation, or transition
width. A very precise specification could ask to achieve the performance goals
with the minimum filter order, or it could call for an arbitrary magnitude
shape, or it might require an FIR filter. Filter design methods differ primarily
in how performance is specified.
To design a filter, the Signal Processing Toolbox software offers two
approaches: object-oriented and non-object oriented. The object-oriented
approach first constructs a filter specification object, fdesign,andthen
invokes an appropriate design method. To illustrate the object-oriented
approach, design and implement a 5–th order lowpass Butterworth filter with
a 3–dB frequency of 200 Hz. Assume a sampling frequency of 1 kHz. Apply
thefiltertoinputdata.
Fs=1000; %Sampling Frequency
time = 0:(1/Fs):1; %time vector
% Data vector
x = cos(2*pi*60*time)+sin(2*pi*120*time)+randn(size(time));
d=fdesign.lowpass('N,F3dB',5,200,Fs); %lowpass filter specification object
% Invoke Butterworth design method
Hd=design(d,'butter');
y=filter(Hd,x);
The non-object oriented approach implements the filter using a function such
as butter and firpm. All of the non-object oriented filter design functions
operate with normalized frequencies. Convert frequency specifications in
Hz to normalized frequency to use these functions. The Signal Processing
Toolbox software defines normalized frequency to be in the closed interval
2-2
Filter Requirements and Specification
[0,1] with 1 denoting πradians/sample. For example, to specify a normalized
frequency of π/2 radians/sample, enter 0.5.
To convert from Hz to normalized frequency, multiply the frequency in Hz by
two and divide by the sampling frequency. To design a 5–th order lowpass
Butterworth filter with a 3–dB frequency of 200 Hz using the non-object
oriented approach, use butter:
Wn = (2*200)/1000; %Convert 3-dB frequency
% to normalized frequency: 0.4*pi rad/sample
[B,A] = butter(5,Wn,'low');
y = filter(B,A,x);
2-3
2Filter Design and Implementation
IIR Filter Design
In this section...
“IIR vs. FIR Filters” on page 2-4
“Classical IIR Filters” on page 2-4
“Other IIR Filters” on page 2-5
“IIR Filter Method Summary” on page 2-5
“Classical IIR Filter Design Using Analog Prototyping” on page 2-6
“Comparison of Classical IIR Filter Types” on page 2-9
IIRvs. FIRFilters
TheprimaryadvantageofIIRfiltersoverFIRfiltersisthattheytypically
meet a given set of specifications with a much lower filter order than a
corresponding FIR filter. Although IIR filters have nonlinear phase, data
processing within MATLAB software is commonly performed “offline,” that
is, the entire data sequence is available prior to filtering. This allows for a
noncausal, zero-phase filtering approach (via the filtfilt function), which
eliminates the nonlinear phase distortion of an IIR filter.
Classical IIR Filters
The classical IIR filters, Butterworth, Chebyshev Types I and II, elliptic, and
Bessel, all approximate the ideal “brick wall” filter in different ways.
This toolbox provides functions to create all these types of classical IIR filters
in both the analog and digital domains (except Bessel, for which only the
analog case is supported), and in lowpass, highpass, bandpass, and bandstop
configurations. For most filter types, you can also find the lowest filter
order that fits a given filter specification in terms of passband and stopband
attenuation, and transition width(s).
2-4
IIR Filter Design
Other IIR Filters
The direct filter design function yulewalk finds a filter with magnitude
response approximating a desired function. This is one way to create a
multiband bandpass filter.
You can also use the parametric modeling or system identification functions
to design IIR filters. These functions are discussed in “Parametric Modeling”
on page 7-13.
The generalized Butterworth design function maxflat is discussed in the
section “Generalized Butterworth Filter Design” on page 2-15.
IIR Filter Method Summary
The following table summarizes the various filter methods in the toolbox and
lists the functionsavailabletoimplementthesemethods.
Toolbox Filters Methods and Available Functions
Filter Method Description Filter Functions
Analog
Prototyping
Using the poles and zeros of a
classical lowpass prototype
filter in the continuous
(Laplace) domain, obtain a
digital filter through frequency
transformation and filter
discretization.
Complete design functions:
besself,butter,cheby1,cheby2,ellip
Order estimation functions:
buttord,cheb1ord,cheb2ord,ellipord
Lowpass analog prototype functions:
besselap,buttap,cheb1ap,cheb2ap,
ellipap
Frequency transformation functions:
lp2bp,lp2bs,lp2hp,lp2lp
Filter discretization functions:
bilinear,impinvar
Direct Design Design digital filter directly
in the discrete time-domain
by approximating a piecewise
linear magnitude response.
yulewalk
2-5
2Filter Design and Implementation
Toolbox Filters Methods and Available Functions (Continued)
Filter Method Description Filter Functions
Generalized
Butterworth
Design
Design lowpass Butterworth
filters with more zeros than
poles.
maxflat
Parametric
Modeling
Find a digital filter that
approximates a prescribed time
or frequency domain response.
(See System Identification
Toolbox™ documentation for
an extensive collection of
parametric modeling tools.)
Time-domain modeling functions:
lpc,prony,stmcb
Frequency-domain modeling functions:
invfreqs,invfreqz
Classical IIR Filter Design Using Analog Prototyping
The principal IIR digital filter design technique this toolbox provides is based
on the conversion of classical lowpass analog filters to their digital equivalents.
The following sections describe how to design filters and summarize the
characteristics of the supported filter types. See “Special Topics in IIR Filter
Design” on page 2-43 for detailed steps on the filter design process.
Complete Classical IIR Filter Design
You can easily create a filter of any order with a lowpass, highpass, bandpass,
or bandstop configuration using the filter design functions.
Filter Design Functions
Filter Type Design Function
Bessel (analog only) [b,a] = besself(n,Wn,options)
[z,p,k] = besself(n,Wn,options)
[A,B,C,D] = besself(n,Wn,options)
Butterworth [b,a] = butter(n,Wn,options)
[z,p,k] = butter(n,Wn,options)
[A,B,C,D] = butter(n,Wn,options)
2-6
IIR Filter Design
Filter Design Functions (Continued)
Filter Type Design Function
Chebyshev Type I [b,a] = cheby1(n,Rp,Wn,options)
[z,p,k] = cheby1(n,Rp,Wn,options)
[A,B,C,D] = cheby1(n,Rp,Wn,options)
Chebyshev Type II [b,a] = cheby2(n,Rs,Wn,options)
[z,p,k] = cheby2(n,Rs,Wn,options)
[A,B,C,D] = cheby2(n,Rs,Wn,options)
Elliptic [b,a] = ellip(n,Rp,Rs,Wn,options)
[z,p,k] = ellip(n,Rp,Rs,Wn,options)
[A,B,C,D] = ellip(n,Rp,Rs,Wn,options)
By default, each of these functions returns a lowpass filter; you need only
specify the desired cutoff frequency Wn in normalized frequency (Nyquist
frequency =1 Hz). For a highpass filter, 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, appending the
string 'stop' for the bandstop configuration.
Herearesomeexampledigitalfilters:
[b,a] = butter(5,0.4); % Lowpass Butterworth
[b,a] = cheby1(4,1,[0.4 0.7]); % Bandpass Chebyshev Type I
[b,a] = cheby2(6,60,0.8,'high'); % Highpass Chebyshev Type II
[b,a] = ellip(3,1,60,[0.4 0.7],'stop'); % Bandstop elliptic
To design an analog filter, perhaps for simulation, use a trailing 's' and
specify cutoff frequencies in rad/s:
[b,a] = butter(5,.4,'s'); % Analog Butterworth filter
All filter design functions return a filter in the transfer function,
zero-pole-gain, or state-space linear system model representation, depending
on how many output arguments are present. In general, you should avoid
2-7
2Filter Design and Implementation
using the transfer function form because numerical problems caused by
roundoff errors can occur. Instead, use the zero-pole-gain form which you can
convert to a second-order section (SOS) form using zp2sos and then use the
SOS form with dfilt to analyze or implement your filter.
Note All classical IIR lowpass filters are ill-conditioned for extremely low
cutoff frequencies. Therefore, instead of designing a lowpass IIR filter with
a very narrow passband, it can be better to design a wider passband and
decimate the input signal.
Designing IIR Filters to Frequency Domain Specifications
This toolbox provides order selection functions that calculate the minimum
filter order that meets a given set of requirements.
Filter Type Order Estimation Function
Butterworth [n,Wn] = buttord(Wp,Ws,Rp,Rs)
Chebyshev Type I [n,Wn] = cheb1ord(Wp, Ws, Rp, Rs)
Chebyshev Type II [n,Wn] = cheb2ord(Wp, Ws, Rp, Rs)
Elliptic [n,Wn] = ellipord(Wp, Ws, Rp, Rs)
These are useful in conjunction with the filter design functions. Suppose you
want a bandpass filter with a passband from 1000 to 2000 Hz, stopbands
starting 500 Hz away on either side, a 10 kHz sampling frequency, at most
1 dB of passband ripple, and at least 60 dB of stopband attenuation. You can
meet these specifications by using the butter function as follows.
[n,Wn] = buttord([1000 2000]/5000,[500 2500]/5000,1,60)
n=
12
Wn =
0.1951 0.4080
[b,a] = butter(n,Wn);
An elliptic filter that meets the same requirements is given by
2-8
IIR Filter Design
[n,Wn] = ellipord([1000 2000]/5000,[500 2500]/5000,1,60)
n=
5
Wn =
0.2000 0.4000
[b,a] = ellip(n,1,60,Wn);
These functions also work with the other standard band configurations, as
well as for analog filters.
Comparison of Classical IIR Filter Types
The toolbox provides five different types of classical IIR filters, each optimal
in some way. This section shows the basic analog prototype form for each and
summarizes major characteristics.
Butterworth Filter
The Butterworth filter provides the best Taylor Series approximation to the
ideal lowpass filter response at analog frequencies Ω=0andΩ=;forany
order N, the magnitude squared response has 2N–1 zero derivatives at these
locations (maximally flat at Ω=0andΩ=). Response is monotonic overall,
decreasing smoothly from Ω=0toΩ=.Hj() /Ω=12
at Ω=1.
2-9
2Filter Design and Implementation
Chebyshev Type I Filter
The Chebyshev Type I filter minimizes the absolute difference between the
ideal and actual frequency response over the entire passband by incorporating
an equal ripple of Rp dB in the passband. Stopband response is maximally
flat. The transition from passband to stopband is more rapid than for the
Butterworth filter. Hj Rp
() /
Ω=
10 20 at Ω=1.
Chebyshev Type II Filter
The Chebyshev Type II filter minimizes the absolute difference between the
ideal and actual frequency response over the entire stopband by incorporating
an equal ripple of Rs dB in the stopband. Passband response is maximally flat.
The stopband does not approach zero as quickly as the type I filter (and
does not approach zero at all for even-valued filter order n). The absence
of ripple in the passband, however, is often an important advantage.
Hj Rs
() /
Ω=
10 20 at Ω=1.
2-10
IIR Filter Design
Elliptic Filter
Elliptic filters are equiripple in both the passband and stopband. They
generally meet filter requirements with the lowest order of any supported
filter type. Given a filter order n, passband ripple Rp in decibels, and
stopband ripple Rs in decibels, elliptic filters minimize transition width.
Hj Rp
() /
Ω=
10 20 at Ω=1.
2-11
2Filter Design and Implementation
Bessel Filter
Analog Bessel lowpass filters have maximally flat group delay at zero
frequency and retain nearly constant group delay across the entire passband.
Filtered signals therefore maintain their waveshapes in the passband
frequency range. Frequency mapped and digital Bessel filters, however, do
not have this maximally flat property; this toolbox supports only the analog
case for the complete Bessel filter design function.
Bessel filters generally require a higherfilterorderthanotherfiltersfor
satisfactory stopband attenuation. Hj() /Ω<12
at Ω= 1 and decreases as
filter order nincreases.
2-12
IIR Filter Design
Note The lowpass filters shown above were created with the analog prototype
functions besselap,buttap,cheb1ap,cheb2ap,andellipap. These functions
find the zeros, poles, and gain of an order nanalog filter of the appropriate
type with cutoff frequency of 1 rad/s. Thecompletefilterdesignfunctions
(besself,butter,cheby1,cheby2,andellip) call the prototyping functions
as a first step in the design process. See “Special Topics in IIR Filter Design”
on page 2-43 for details.
To create similar plots, use n=5and, as needed, Rp =0.5 and Rs =20.For
example, to create the elliptic filter plot:
[z,p,k] = ellipap(5,0.5,20);
w = logspace(-1,1,1000);
h = freqs(k*poly(z),poly(p),w);
semilogx(w,abs(h)), grid
Direct IIR Filter Design
This toolbox uses the term direct methods to describe techniques for IIR
design that find a filter based on specifications in the discrete domain. Unlike
the analog prototyping method, direct design methods are not constrained
to the standard lowpass, highpass, bandpass, or bandstop configurations.
Rather, these functions design filters with an arbitrary, perhaps multiband,
frequency response. This section discusses the yulewalk function, which
is intended specifically for filter design; “Parametric Modeling” on page
7-13 discusses other methods that may also be considered direct, such as
Prony’s method, Linear Prediction, the Steiglitz-McBride method, and inverse
frequency design.
The yulewalk function designs recursive IIR digital filters by fitting a
specified frequency response. yulewalk’s name reflects its method for
finding the filter’s denominator coefficients: it finds the inverse FFT of
the ideal desired magnitude-squared response and solves the modified
Yule-Walker equations using the resulting autocorrelation function samples.
The statement
[b,a] = yulewalk(n,f,m)
2-13
2Filter Design and Implementation
returns row vectors band acontaining the n+1 numerator and denominator
coefficients of the order nIIR filter whose frequency-magnitude characteristics
approximate those given in vectors fand m.fis a vector of frequency points
ranging from 0 to 1, where 1 represents the Nyquist frequency. mis a vector
containing the desired magnitude response at the points in f.fand m
can describe any piecewise linear shape magnitude response, including a
multiband response. The FIR counterpart of this function is fir2,whichalso
designs a filter based on an arbitrary piecewise linear magnitude response.
See “FIR Filter Design” on page 2-17 for details.
Note that yulewalk does not accept phase information, and no statements are
made about the optimality of the resulting filter.
Design a multiband filter with yulewalk, and plot the desired and actual
frequency response:
m=[0011001100];
f = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1];
[b,a] = yulewalk(10,f,m);
[h,w] = freqz(b,a,128)
plot(f,m,w/pi,abs(h))
2-14
IIR Filter Design
Generalized Butterworth Filter Design
The toolbox function maxflat enables you to design generalized Butterworth
filters, that is, Butterworth filters with differing numbers of zeros and poles.
This is desirable in some implementations where poles are more expensive
computationally than zeros. maxflat is just like the butter function, except
that it you can specify two orders (one for the numerator and one for the
denominator) instead of just one. These filters are maximally flat.This
means that the resulting filter is optimal for any numerator and denominator
orders, with the maximum number of derivatives at 0 and the Nyquist
frequency ω=πboth set to 0.
For example, when the two orders are the same, maxflat is the same as
butter:
[b,a] = maxflat(3,3,0.25)
b=
0.0317 0.0951 0.0951 0.0317
a=
1.0000 -1.4590 0.9104 -0.1978
[b,a] = butter(3,0.25)
b=
0.0317 0.0951 0.0951 0.0317
a=
1.0000 -1.4590 0.9104 -0.1978
However, maxflat is more versatile because it allows you to design a filter
with more zeros than poles:
[b,a] = maxflat(3,1,0.25)
b=
0.0950 0.2849 0.2849 0.0950
a=
1.0000 -0.2402
The third input to maxflat is the half-power frequency, a frequency between
0 and 1 with a desired magnitude response of 12/.
You can also design linear phase filters that have the maximally flat property
using the 'sym' option:
2-15
2Filter Design and Implementation
maxflat(4,'sym',0.3)
ans =
0.0331 0.2500 0.4337 0.2500 0.0331
For complete details of the maxflat algorithm, see Selesnick and Burrus [2].
2-16
FIR Filter Design
FIR Filter Design
In this section...
“FIR vs. IIR Filters” on page 2-17
“FIR Filter Summary” on page 2-18
“Linear Phase Filters” on page 2-18
“Windowing Method” on page 2-20
“Multiband FIR Filter Design with Transition Bands” on page 2-24
“Constrained Least Squares FIR Filter Design” on page 2-31
“Arbitrary-Response Filter Design” on page 2-37
FIRvs. IIRFilters
Digital filters with finite-duration impulse response (all-zero, or FIR filters)
have both advantages and disadvantages compared to infinite-duration
impulseresponse(IIR)filters.
FIR filters have the following primary advantages:
They can have exactly linear phase.
They are always stable.
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.
2-17
2Filter Design and Implementation
FIR Filter Summary
FIR Filters
Filter Design
Method Description Filter Functions
Windowing Apply window to truncated inverse
Fourier transform of desired “brick
wall” filter
fir1,fir2,
kaiserord
Multiband
with Transition
Bands
Equiripple or least squares approach
over sub-bands of the frequency range
firls,firpm,
firpmord
Constrained
Least Squares
Minimize squared integral error over
entire frequency range subject to
maximum error constraints
fircls,fircls1
Arbitrary
Response
Arbitrary responses, including
nonlinear phase and complex filters
cfirpm
Raised Cosine Lowpass response with smooth,
sinusoidal transition
firrcos
Linear Phase Filters
Except for cfirpm, all of the FIR filter design functions design linear phase
filters only. The filter coefficients, or “taps,” of such filters obey either an even
or odd symmetry relation. Depending on this symmetry, and on whether the
ordernof the filter is even or odd, a linear phase filter (stored in length n+1
vector b) has certain inherent restrictions on its frequency response.
2-18
FIR Filter Design
Linear
Phase
Filter Type
Filter
Order Symmetry of Coefficients
Response
H(f), f =0
Response
H(f), f =1
(Nyquist)
Type I Even even:
bk bn k k n( ) ( ), ,...,=+− = +211
No
restriction
No
restriction
Type II Odd even:
bk bn k k n( ) ( ), ,...,=+− = +211
No
restriction
H(1) =0
Type III Even odd:
bk bn k k n( ) ( ), ,...,=− + − = +211
H(0) =0H(1) =0
Type IV Odd odd:
bk bn k k n( ) ( ), ,...,=− + − = +211
H(0) =0No
restriction
The phase delay and group delay of linear phase FIR filters are equal and
constant over the frequency band. For an order nlinear phase FIR filter, the
group delay is n/2, and the filtered signal is simply delayed by n/2 time steps
(and the magnitude of its Fourier transform is scaled by the filter’s magnitude
response). This property preserves the wave shape of signals in the passband;
that is, there is no phase distortion.
The functions fir1,fir2,firls,firpm,fircls,fircls1,andfirrcos 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. cfirpm can design any type of linear phase filter,
and nonlinear phase filters as well.
Note 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 nin these cases, fir1 adds 1 to the order
and returns a type I filter.
2-19
2Filter Design and Implementation
Windowing Method
Consider the ideal, or “brick wall,” digital lowpass filter with a cutoff
frequency of ω0rad/s. This filter has magnitude 1 at all frequencies with
magnitude less than ω0, and magnitude 0 at frequencies with magnitude
between ω0and π. Its impulse response sequence h(n)is
hn H e d e d n
n
jn jn
() ( ) sin( )
===
−−
∫∫
1
2
1
20
00
πωω
πωω
π
π
πωω
ω
ω
This filter is not implementable since its impulse response is infinite and
noncausal. To create a finite-duration impulse response, truncate it by
applying a window. By retaining the central section of impulse response in
this truncation, you obtain a linear phase FIR filter. For example, a length 51
filter with a lowpass cutoff frequency ω0of 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 thatbestapproximatestheideallowpass
filter, in the integrated least squares sense. The following command displays
the filter’s frequency response in FVTool:
fvtool(b,1)
2-20
FIR Filter Design
Note that the y-axis shown in the figure below is in Magnitude Squared.
You can set this by right-clicking on the axis label and selecting Magnitude
Squared from the menu.
Ringing and ripples occur in the response, especially near the band edge.
This “Gibbs effect” 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 smoothing in the frequency domain.
Apply a length 51 Hamming window to the filter and display the result using
FVTool:
b = 0.4*sinc(0.4*(-25:25));
b = b.*hamming(51)';
fvtool(b,1)
2-21
2Filter Design and Implementation
Note that the y-axis shown in the figure below is in Magnitude Squared.
You can set this by right-clicking on the axis label and selecting Magnitude
Squared from the menu.
Using a Hamming window greatly reduces the ringing. This improvement is
at the expense of transition width (the windowed version takes longer to ramp
from passband to stopband) and optimality (the windowed version does not
minimize the integrated squared error).
The functions fir1 and fir2 are based on this windowing process. Given a
filter order and description of an ideal desired filter, these functions return a
windowed inverse Fourier transform of that ideal filter. Both use a Hamming
window by default, but they accept any window function. See “Windows” on
page 7-2 for an overview of windows and their properties.
2-22
FIR Filter Design
Standard Band FIR Filter Design: fir1
fir1 implements the classical method of windowed linear phase FIR
digital filter design. It resembles the IIR filter design functions in that it
is formulated to design filters in standard band configurations: lowpass,
bandpass, highpass, and bandstop.
The statements
n = 50;
Wn = 0.4;
b = fir1(n,Wn);
create row vector bcontaining the coefficients of the order n
Hamming-windowed filter. This is a lowpass, linear phase FIR filter with
cutoff frequency Wn.Wn is a number between 0 and 1, where 1 corresponds to
the Nyquist frequency, half the sampling frequency. (Unlike other methods,
here Wn corresponds to the 6 dB point.) 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.
Kaiser Window Order Estimation. The kaiserord function estimates the
filter order, cutoff frequency, and Kaiser window beta parameter needed to
meet a given set of specifications. Given a vector of frequency band edges and
a corresponding vector of magnitudes, as well as maximum allowable ripple,
kaiserord returns appropriate input parameters for the fir1 function.
Multiband FIR Filter Design: fir2
The fir2 function also designs windowed FIR filters, but with an arbitrarily
shaped piecewise linear frequency response. This is in contrast to fir1,which
only designs filters in standard lowpass, highpass, bandpass, and bandstop
configurations.
The commands
2-23
2Filter Design and Implementation
n = 50;
f=[0.4.51];
m=[1 1 00];
b = fir2(n,f,m);
return row vector bcontaining the n+1 coefficients of the order nFIR filter
whose frequency-magnitude characteristics match those given by vectors f
and m.fis a vector of frequency points ranging from 0 to 1, where 1 represents
the Nyquist frequency. mis a vector containing the desired magnitude
response at the points specified in f. (The IIR counterpart of this function
is yulewalk, which also designs filters based on arbitrary piecewise linear
magnitude responses. See “IIR Filter Design” on page 2-4 for details.)
Multiband FIR Filter Design with Transition Bands
The firls and firpm functions provide a more general means of specifying
the ideal desired filter than the fir1 and fir2 functions. These functions
design Hilbert transformers, differentiators, and other filters with odd
symmetric coefficients (type III and type IV linear phase). They also let you
include transition or “don’t care” regions in which the error is not minimized,
and perform band dependent weighting of the minimization.
The firls function is an extension of the fir1 and fir2 functions in that
it minimizes the integral of the square of the error between the desired
frequency response and the actual frequency response.
The firpm function implements the Parks-McClellan algorithm, which uses
the Remez exchange algorithm and Chebyshev approximation theory to design
filters with optimal fits between the desired and actual frequency responses.
The filters are optimal in the sense that they minimize the maximum error
between the desired frequency response and the actual frequency response;
they are sometimes called minimax filters. Filters designed in this way
exhibit an equiripple behavior in their frequency response, and hence are also
known as equiripple filters. The Parks-McClellan FIR filter design algorithm
is perhaps the most popular and widely used FIR filter design methodology.
The syntax for firls and firpm isthesame;theonlydifferenceistheir
minimization schemes. The next example shows how filters designed with
firls and firpm reflect these different schemes.
2-24
FIR Filter Design
Basic Configurations
The default mode of operation of firls and firpm is to design type I or type
II linear phase filters, depending on whether the order you desire is even or
odd, respectively. A lowpass example with approximate amplitude 1 from 0 to
0.4 Hz, and approximate amplitude 0 from 0.5 to 1.0 Hz is
n = 20; % Filter order
f = [0 0.4 0.5 1]; % Frequency band edges
a = [1 1 0 0]; % Desired amplitudes
b = firpm(n,f,a);
From 0.4 to 0.5 Hz, firpm performs no error minimization; this is a transition
band or “don’t care” region. A transition band minimizes the error more in
thebandsthatyoudocareabout,attheexpenseofaslowertransitionrate.
In this way, these types of filters have an inherent trade-off similar to FIR
design by windowing.
Tocompareleastsquarestoequiripplefilterdesign,usefirls to create
a similar filter. Type
bb = firls(n,f,a);
and compare their frequency responses using FVTool:
fvtool(b,1,bb,1)
2-25
2Filter Design and Implementation
Note that the y-axis shown in the figure below is in Magnitude Squared.
You can set this by right-clicking on the axis label and selecting Magnitude
Squared from the menu.
The filter designed with firpm exhibits equiripple behavior. Also note that
the firls filter has a better response over most of the passband and stopband,
but at the band edges (f=0.4 and f=0.5),theresponseisfurtherawayfrom
the ideal than the firpm filter. This shows that the firpm filter’s maximum
error over the passband and stopband is smaller and, in fact, it is the smallest
possible for this band edge configuration and filter length.
Think of frequency bands as lines over short frequency intervals. firpm and
firls use this scheme to represent any piecewise linear desired function with
any transition bands. firls and firpm design lowpass, highpass, bandpass,
and bandstop filters; a bandpass example is
f = [0 0.3 0.4 0.7 0.8 1]; % Band edges in pairs
2-26
FIR Filter Design
a = [0 0 1 1 0 0]; % Bandpass filter amplitude
Technically, these fand avectors define five bands:
Two stopbands, from 0.0 to 0.3 and from 0.8 to 1.0
A passband from 0.4 to 0.7
Two transition bands, from 0.3 to 0.4 and from 0.7 to 0.8
Example highpass and bandstop filters are
f = [0 0.7 0.8 1]; % Band edges in pairs
a = [0 0 1 1]; % Highpass filter amplitude
f = [0 0.3 0.4 0.5 0.8 1]; % Band edges in pairs
a = [1 1 0 0 1 1]; % Bandstop filter amplitude
An example multiband bandpass filter is
f = [0 0.1 0.15 0.25 0.3 0.4 0.45 0.55 0.6 0.7 0.75 0.85 0.9 1];
a=[11001100110011];
Another possibility is a filter that has as a transition region the line
connecting the passband with the stopband; this can help control “runaway”
magnitude response in wide transition regions:
f = [0 0.4 0.42 0.48 0.5 1];
a = [1 1 0.8 0.2 0 0]; % Passband, linear transition,
% stopband
The Weight Vector
Both firls and firpm allow you to place more or less emphasis on minimizing
the error in certain frequency bands relative to others. To do this, specify a
weight vector following the frequencyandamplitudevectors. Anexample
lowpass equiripple filter with 10 times less ripple in the stopband than the
passband is
n = 20; % Filter order
f = [0 0.4 0.5 1]; % Frequency band edges
a = [1 1 0 0]; % Desired amplitudes
2-27
2Filter Design and Implementation
w = [1 10]; % Weight vector
b = firpm(n,f,a,w);
A legal weight vector is always half the length of the fand avectors; there
must be exactly one weight per band.
Anti-Symmetric Filters / Hilbert Transformers
When called with a trailing 'h' or 'Hilbert' option, firpm and firls design
FIR filters with odd symmetry, that is, type III (for even order) or type IV
(for odd order) linear phase filters. An ideal Hilbert transformer has this
anti-symmetry property and an amplitude of 1 across the entire frequency
range. Try the following approximate Hilbert transformers and plot them
using FVTool:
b = firpm(21,[0.05 1],[1 1],'h'); % Highpass Hilbert
bb = firpm(20,[0.05 0.95],[1 1],'h'); % Bandpass Hilbert
fvtool(b,1,bb,1)
2-28
FIR Filter Design
You can find the delayed Hilbert transform of a signal xby passing it through
these filters.
fs = 1000; % Sampling frequency
t = (0:1/fs:2)'; % Two second time vector
x = sin(2*pi*300*t); % 300 Hz sine wave example signal
xh = filter(bb,1,x); % Hilbert transform of x
The analytic signal corresponding to xis the complex signal that has xas its
real part and the Hilbert transform of xas its imaginary part. For this FIR
method (an alternative to the hilbert function), you must delay xby half the
filter order to create the analytic signal:
xd = [zeros(10,1); x(1:length(x)-10)]; % Delay 10 samples
xa = xd + j*xh; % Analytic signal
2-29
2Filter Design and Implementation
This method does not work directly for filters of odd order, which require a
noninteger delay. In this case, the hilbert function, described in “Specialized
Transforms” on page 7-36, estimates the analytic signal. Alternatively, use
the resample function to delay the signal by a noninteger number of samples.
Differentiators
Differentiation of a signal in the time domain is equivalent to multiplication
of the signal’s Fourier transform by an imaginary ramp function. That is, to
differentiate a signal, pass it through a filter that has a response H(ω)=jω.
Approximate the ideal differentiator (with a delay) using firpm or firls with
a'd' or 'differentiator' option:
b = firpm(21,[0 1],[0 pi],'d');
For a type III filter, the differentiation band should stop short of the Nyquist
frequency, and the amplitude vector must reflect that change to ensure the
correct slope:
bb = firpm(20,[0 0.9],[0 0.9*pi],'d');
In the 'd' mode, firpm weights the error by 1/ωin nonzero amplitude bands
to minimize the maximum relative error. firls weights the error by (1/ω)2in
nonzero amplitude bands in the 'd' mode.
2-30
FIR Filter Design
The following plots show the magnitude responses for the differentiators
above.
fvtool(b,1,bb,1)
Constrained Least Squares FIR Filter Design
The Constrained Least Squares (CLS) FIR filter design functions implement
a technique that enables you to design FIR filters without explicitly defining
the transition bands for the magnitude response. The ability to omit the
specification of transition bands is useful in several situations. For example,
it may not be clear where a rigidly defined transition band should appear if
noise and signal information appear together in the same frequency band.
Similarly, it may make sense to omit the specification of transition bands if
they appear only to control the results of Gibbs phenomena that appear in
the filter’s response. See Selesnick, Lang, and Burrus [2] for discussion of
this method.
2-31
2Filter Design and Implementation
Instead of defining passbands, stopbands, and transition regions, the CLS
method accepts a cutoff frequency (for the highpass, lowpass, bandpass, or
bandstop cases), or passband and stopband edges (for multiband cases), for
the desired response. In this way, the CLS method defines transition regions
implicitly, rather than explicitly.
The key feature of the CLS method is that it enables you to define upper
and lower thresholds that contain the maximum allowable ripple in the
magnitude response. Given this constraint, the technique applies the least
square error minimization technique over the frequency range of the filter’s
response, instead of over specific bands. The error minimization includes
any areas of discontinuity in the ideal, “brick wall” response. An additional
benefit is that the technique enables you to specify arbitrarily small peaks
resulting from Gibbs’ phenomena.
There are two toolbox functions that implement this design technique.
Description Function
Constrained least square multiband FIR filter design fircls
Constrained least square filter design for lowpass and
highpass linear phase filters
fircls1
For details on the calling syntax for these functions, see their reference
descriptions in the Function Reference.
Basic Lowpass and Highpass CLS Filter Design
The most basic of the CLS design functions, fircls1, uses this technique to
design lowpass and highpass FIR filters. As an example, consider designing a
filter with order 61 impulse response and cutoff frequency of 0.3 (normalized).
Further, define the upper and lower bounds that constrain the design process
as:
Maximum passband deviation from 1 (passband ripple) of 0.02.
Maximum stopband deviation from 0 (stopband ripple) of 0.008.
2-32
FIR Filter Design
To approach this design problem using fircls1, use the following commands:
n = 61;
wo = 0.3;
dp = 0.02;
ds = 0.008;
h = fircls1(n,wo,dp,ds);
fvtool(h,1)
Note that the y-axisshownbelowisinMagnitude Squared. You can set
this by right-clicking on the axis label and selecting Magnitude Squared
from the menu.
2-33
2Filter Design and Implementation
Multiband CLS Filter Design
fircls uses the same technique to design FIR filters with a desired piecewise
constant magnitude response. In this case, you can specify a vector of band
edges and a corresponding vector of band amplitudes. In addition, you can
specify the maximum amount of ripple for each band.
For example, assume the specifications for a filter call for:
From 0 to 0.3 (normalized): amplitude 0, upper bound 0.005, lower
bound –0.005
From 0.3 to 0.5: amplitude 0.5, upper bound 0.51, lower bound 0.49
From 0.5 to 0.7: amplitude 0, upper bound 0.03, lower bound –0.03
From 0.7 to 0.9: amplitude 1, upper bound 1.02, lower bound 0.98
From 0.9 to 1: amplitude 0, upper bound 0.05, lower bound –0.05
Design a CLS filter with impulse response order 129 that meets these
specifications:
n = 129;
f = [0 0.3 0.5 0.7 0.9 1];
a = [0 0.5 0 1 0];
up = [0.005 0.51 0.03 1.02 0.05];
lo = [-0.005 0.49 -0.03 0.98 -0.05];
h = fircls(n,f,a,up,lo);
fvtool(h,1)
2-34
FIR Filter Design
Note that the y-axisshownbelowisinMagnitude Squared. You can set
this by right-clicking on the axis label and selecting Magnitude Squared
from the menu.
Weighted CLS Filter Design
Weighted CLS filter design lets you design lowpass or highpass FIR filters
with relative weighting of the error minimization in each band. The fircls1
function enables you to specify the passband and stopband edges for the least
squares weighting function, as well as a constant kthat specifies the ratio of
the stopband to passband weighting.
Forexample,considerspecificationsthatcallforanFIRfilterwithimpulse
response order of 55 and cutoff frequency of 0.3 (normalized). Also assume
maximum allowable passband ripple of 0.02 and maximum allowable
stopband ripple of 0.004. In addition, add weighting requirements:
2-35
2Filter Design and Implementation
Passband edge for the weight function of 0.28 (normalized)
Stopband edge for the weight function of 0.32
Weight error minimization 10 times as much in the stopband as in the
passband
To approach this using fircls1,type
n = 55;
wo = 0.3;
dp = 0.02;
ds = 0.004;
wp = 0.28;
ws = 0.32;
k = 10;
h = fircls1(n,wo,dp,ds,wp,ws,k);
fvtool(h,1)
2-36
FIR Filter Design
Note that the y-axisshownbelowisinMagnitude Squared. You can set
this by right-clicking on the axis label and selecting Magnitude Squared
from the menu.
Arbitrary-Response Filter Design
The cfirpm filter design function provides a tool for designing FIR filters
with arbitrary complex responses. It differs from the other filter design
functions in how the frequency response of the filter is specified: it accepts the
name of a function which returns the filter response calculated over a grid of
frequencies. This capability makes cfirpm a highly versatile and powerful
technique for filter design.
This design technique may be used to produce nonlinear-phase FIR filters,
asymmetric frequency-response filters (with complex coefficients), or more
symmetric filters with custom frequency responses.
2-37
2Filter Design and Implementation
The design algorithm optimizes the Chebyshev (or minimax) error using
an extended Remez-exchange algorithm for an initial estimate. If this
exchange method fails to obtain the optimal filter, the algorithm switches
to an ascent-descent algorithm that takes over to finish the convergence to
the optimal solution.
Multiband Filter Design
Consider a multiband filter with the following special frequency-domain
characteristics.
Band Amplitude
Optimization
Weighting
[–1 –0.5] [5 1] 1
[–0.4 +0.3] [2 2] 10
[+0.4 +0.8] [2 1] 5
A linear-phase multiband filter may be designed using the predefined
frequency-response function multiband,asfollows:
b = cfirpm(38, [-1 -0.5 -0.4 0.3 0.4 0.8], ...
{'multiband', [5 1 2 2 2 1]}, [1 10 5]);
For the specific case of a multiband filter, we can use a shorthand filter design
notation similar to the syntax for firpm:
b = cfirpm(38,[-1 -0.5 -0.4 0.3 0.4 0.8], ...
[5 1 2 2 2 1], [1 10 5]);
As with firpm, a vector of band edges is passed to cfirpm. This vector defines
the frequency bands over which optimization is performed; note that there are
two transition bands, from –0.5 to –0.4 and from 0.3 to 0.4.
In either case, the frequency response is obtained and plotted using linear
scale in FVTool:
fvtool(b,1)
2-38
FIR Filter Design
Note that the range of data shown below is (-Fs/2,Fs/2). You can set this
range by changing the x-axis units to Frequency (Fs = 1 Hz).
2-39
2Filter Design and Implementation
The filter response for this multiband filter is complex, which is expected
because of the asymmetry in the frequency domain. The impulse response,
which you can select from the FVTool toolbar, is shown below.
Filter Design with Reduced Delay
Consider the design of a 62-tap lowpass filter with a half-Nyquist cutoff. If
we specify a negative offset value to the lowpass filter design function, the
group delay offset for the design is significantly less than that obtained for a
standard linear-phase design. This filter design may be computed as follows:
b = cfirpm(61,[0 0.5 0.55 1],{'lowpass',-16});
The resulting magnitude response is
fvtool(b,1)
2-40
FIR Filter Design
Note that the range of data in this plot is (-Fs/2,Fs/2),whichyoucanset
changing the x-axis units to Frequency.They-axisisinMagnitudeSquared,
which you can set by right-clicking on the axis label and selecting Magnitude
Squared from the menu.
2-41
2Filter Design and Implementation
The group delay of the filter reveals that the offset has been reduced from
N/2 to N/2-16 (i.e., from 30.5 to 14.5). Now, however, the group delay is no
longer flat in the passband region. To create this plot, click the Group Delay
button on the toolbar.
If we compare this nonlinear-phase filter to a linear-phase filter that has
exactly 14.5 samples of group delay, the resulting filter is of order 2*14.5, or
29. Using b = cfirpm(29,[0 0.5 0.55 1],'lowpass'), the passband and
stopband ripple is much greater for the order 29 filter. These comparisons
can assist you in deciding which filter is more appropriate for a specific
application.
2-42
Special Topics in IIR Filter Design
Special Topics in IIR Filter Design
In this section...
“Classic IIR Filter Design” on page 2-43
“Analog Prototype Design” on page 2-44
“Frequency Transformation” on page 2-44
“Filter Discretization” on page 2-46
Classic IIR Filter Design
The classic IIR filter design technique includes the following steps.
1Find an analog lowpass filter with cutoff frequency of 1 and translate this
prototype filter to the desired band configuration
2Transform the filter to the digital domain.
3Discretize the filter.
The toolbox provides functions for each of these steps.
Design Task Available functions
Analog lowpass
prototype
buttap,cheb1ap,besselap,ellipap,cheb2ap
Frequency
transformation
lp2lp,lp2hp,lp2bp,lp2bs
Discretization bilinear,impinvar
Alternatively, the butter,cheby1,cheb2ord,ellip,andbesself functions
perform all steps of the filter design and the buttord,cheb1ord,cheb2ord,
and ellipord functions provide minimum order computation for IIR filters.
These functions are sufficient for many design problems, and the lower level
functions are generally not needed. But if you do have an application where
you need to transform the band edges of an analog filter, or discretize a
rational transfer function, this section describes the tools with which to do so.
2-43
2Filter Design and Implementation
Analog PrototypeDesign
This toolbox provides a number of functions to create lowpass analog
prototype filters with cutoff frequency of 1, the first step in the classical
approach to IIR filter design.
The table below summarizes the analog prototype design functions for each
supported filter type; plots for each type are shown in “IIR Filter Design”
on page 2-4.
Filter Type Analog Prototype Function
Bessel [z,p,k] = besselap(n)
Butterworth [z,p,k] = buttap(n)
Chebyshev Type I [z,p,k] = cheb1ap(n,Rp)
Chebyshev Type II [z,p,k] = cheb2ap(n,Rs)
Elliptic [z,p,k] = ellipap(n,Rp,Rs)
Frequency Transformation
The second step in the analog prototyping design technique is the frequency
transformation of a lowpass prototype. The toolbox provides a set of functions
to transform analog lowpass prototypes (with cutoff frequency of 1 rad/s)
into bandpass, highpass, bandstop, and lowpass filters of the desired cutoff
frequency.
Frequency
Transformation Transformation Function
Lowpass to lowpass
=ss/
ω
0
[numt,dent] = lp2lp (num,den,Wo)
[At,Bt,Ct,Dt] =lp2lp (A,B,C,D,Wo)
Lowpass to highpass
=ss
ω
0
[numt,dent] = lp2hp (num,den,Wo)
[At,Bt,Ct,Dt] =lp2hp (A,B,C,D,Wo)
2-44
Special Topics in IIR Filter Design
Frequency
Transformation Transformation Function
Lowpass to bandpass
=+
sB
s
s
ωω
ω
ω
00
2
0
1(/ )
/
[numt,dent] = lp2bp (num,den,Wo,Bw)
[At,Bt,Ct,Dt] =lp2bp (A,B,C,D,Wo,Bw)
Lowpass to bandstop
=+
sBs
s
ω
ω
ω
ω
0
0
021
/
(/ )
[numt,dent] = lp2bs (num,den,Wo,Bw)
[At,Bt,Ct,Dt] =lp2bs( A,B,C,D,Wo,Bw)
As shown, all of the frequency transformation functions can accept two linear
system models: transfer function and state-space form. For the bandpass
and bandstop cases
ωωω
012
=
and
B
ω
ωω
=−
21
where ω1is the lower band edge and ω2is the upper band edge.
The frequency transformation functions perform frequency variable
substitution. In the case of lp2bp and lp2bs, this is a second-order
substitution, so the output filter is twice the order of the input. For lp2lp and
lp2hp, the output filter is the same order as the input.
To begin designing an order 10 bandpass Chebyshev Type I filter with a value
of 3 dB for passband ripple, enter
[z,p,k] = cheb1ap(10,3);
Outputs z,p,andkcontain the zeros, poles, and gain of a lowpass analog
filter with cutoff frequency Ωcequal to 1 rad/s. Use the lp2bp function to
transform this lowpass prototype to a bandpass analog filter with band edges
2-45
2Filter Design and Implementation
Ω1=π/5 and Ω2=π. First, convert the filter to state-space form so the lp2bp
function can accept it:
[A,B,C,D] = zp2ss(z,p,k); % Convert to state-space form.
Now, find the bandwidth and center frequency, and call lp2bp:
u1 = 0.1*2*pi; u2 = 0.5*2*pi; % In radians per second
Bw = u2-u1;
Wo = sqrt(u1*u2);
[At,Bt,Ct,Dt] = lp2bp(A,B,C,D,Wo,Bw);
Finally, calculate the frequency response and plot its magnitude:
[b,a] = ss2tf(At,Bt,Ct,Dt); % Convert to TF form
w = linspace(0.01,1,500)*2*pi; % Generate frequency vector
h = freqs(b,a,w); % Compute frequency response
semilogy(w/2/pi,abs(h)), grid % Plot log magnitude vs. freq
xlabel('Frequency (Hz)');
Filter Discretization
Thethirdstepintheanalogprototyping technique is the transformation of
thefiltertothediscrete-timedomain. The toolbox provides two methods for
this: the impulse invariant and bilinear transformations. The filter design
2-46
Special Topics in IIR Filter Design
functions butter,cheby1,cheby2,andellip use the bilinear transformation
for discretization in this step.
Analog
to Digital
Transformation Transformation Function
Impulse
invariance
[numd,dend] = impinvar (num,den,fs)
Bilinear
transform
[zd,pd,kd] = bilinear (z,p,k,fs,Fp)
[numd,dend] = bilinear (num,den,fs,Fp)
[Ad,Bd,Cd,Dd] = bilinear (At,Bt,Ct,Dt,fs,Fp)
Impulse Invariance
The toolbox function impinvar creates a digital filter whose impulse response
is the samples of the continuous impulse response of an analog filter. This
function works only on filters in transfer function form. For best results, the
analog filter should have negligible frequency content above half the sampling
frequency, because such high frequency content is aliased into lower bands
upon sampling. Impulse invariance works for some lowpass and bandpass
filters, but is not appropriate for highpass and bandstop filters.
Design a Chebyshev Type I filter and plot its frequency and phase response
using FVTool:
[bz,az] = impinvar(b,a,2);
fvtool(bz,az)
Click the Magnitude and Phase Response toolbar button.
2-47
2Filter Design and Implementation
Impulse invariance retains the cutoff frequencies of 0.1 Hz and 0.5 Hz.
Bilinear Transformation
The bilinear transformation is a nonlinear mapping of the continuous domain
to the discrete domain; it maps the s-plane into the z-plane by
Hz Hs sk
z
z
() ()==
+
1
1
Bilinear transformation maps the jΩ-axis of the continuous domain to the
unit circle of the discrete domain according to
ω
=
21
tan Ω
k
2-48
Special Topics in IIR Filter Design
The toolbox function bilinear implements this operation, where the
frequency warping constant kis equal to twice the sampling frequency (2*fs)
by default, and equal to 2ππfff
pps
tan
()
if you give bilinear a trailing
argument that represents a “match” frequency Fp. If a match frequency Fp
(in hertz) is present, bilinear maps the frequency Ω=2πfp(in rad/s) to the
same frequency in the discrete domain, normalized to the sampling rate:
ω=2πfp/fs(in rad/sample).
The bilinear function can perform this transformation on three different
linear system representations: zero-pole-gain, transfer function, and
state-space form. Try calling bilinear with the state-space matrices that
describe the Chebyshev Type I filter from the previous section, using a
sampling frequency of 2 Hz, and retaining the lower band edge of 0.1 Hz:
[Ad,Bd,Cd,Dd] = bilinear(At,Bt,Ct,Dt,2,0.1);
The frequency response of the resulting digital filter is
[bz,az] = ss2tf(Ad,Bd,Cd,Dd); % Convert to TF
fvtool(bz,az)
2-49
2Filter Design and Implementation
Click the Magnitude and Phase Response toolbar button.
The lower band edge is at 0.1 Hz as expected. Notice, however, that the
upper band edge is slightly less than 0.5 Hz, although in the analog domain
it was exactly 0.5 Hz. This illustrates the nonlinear nature of the bilinear
transformation. To counteract this nonlinearity, it is necessary to create
analog domain filters with “prewarped” band edges, which map to the correct
locations upon bilinear transformation. Here the prewarped frequencies u1
and u2 generate Bw and Wo for the lp2bp function:
fs = 2; % Sampling frequency (hertz)
u1 = 2*fs*tan(0.1*(2*pi/fs)/2); % Lower band edge (rad/s)
u2 = 2*fs*tan(0.5*(2*pi/fs)/2); % Upper band edge (rad/s)
Bw = u2 - u1; % Bandwidth
Wo = sqrt(u1*u2); % Center frequency
[At,Bt,Ct,Dt] = lp2bp(A,B,C,D,Wo,Bw);
2-50
Special Topics in IIR Filter Design
A digital bandpass filter with correct band edges 0.1 and 0.5 times the
Nyquist frequency is
[Ad,Bd,Cd,Dd] = bilinear(At,Bt,Ct,Dt,fs);
The example bandpass filters from the last two sections could also be created
in one statement using the complete IIR design function cheby1.Forinstance,
an analog version of the example Chebyshev filter is
[b,a] = cheby1(5,3,[0.1 0.5]*2*pi,'s');
Note that the band edges are in rad/s for analog filters, whereas for the digital
case, frequency is normalized:
[bz,az] = cheby1(5,3,[0.1 0.5]);
All of the complete design functions call bilinear internally. They prewarp
the band edges as needed to obtain the correct digital filter.
2-51
2Filter Design and Implementation
Filtering Data With Signal Processing Toolbox Software
In this section...
“Lowpass FIR Filter — Window Method” on page 2-52
“Bandpass Filters — Minimum-Order FIR and IIR Systems” on page 2-56
“Zero-Phase Filtering” on page 2-65
Lowpass FIR Filter — Window Method
These examples show you how to design and implement an FIR filter using
the command line functions: fir1 and fdesign.lowpass, and the interactive
tool fdatool.
To filter the input signal, these examples use the filter command. The
examples in “Zero-Phase Filtering” on page 2-65 show you how to implement
zero-phase filtering with filtfilt.
Createasignaltouseintheexamples. Thesignalisa100-Hzsinewavein
additive N(0,1/4) white Gaussian noise. Set the random number generator to
the default state for reproducible results.
rng default;
Fs = 1000;
t = linspace(0,1,Fs);
x = cos(2*pi*100*t)+0.5*randn(size(t));
The filter design is an FIR lowpass filter with order equal to 20 and a cutoff
frequency of 150 Hz. Use a Kasier window with length one sample greater
than the filter order and β=3. See kaiser fordetailsontheKaiserwindow.
Use fir1 to design the filter. fir1 requires normalized frequencies in the
interval [0,1], where 1 is (1)πradians/sample. To use fir1,youmustconvert
all frequency specifications to normalized frequencies.
Design the filter and view the filter’s magnitude response.
fc = 150;
Wn = (2/Fs)*fc;
2-52
Filtering Data With Signal Processing Toolbox™ Software
b = fir1(20,Wn,'low',kaiser(21,3));
fvtool(b,1,'Fs',Fs);
Apply the filter to the signal and plot the result for the first ten periods of
the 100-Hz sinusoid.
y = filter(b,1,x);
plot(t(1:100),x(1:100),'k');
hold on;
plot(t(1:100),y(1:100),'r','linewidth',2);
legend('Original Signal','Filtered Data','Location','SouthEast');
xlabel('Seconds'); ylabel('Amplitude');
Designthesamefilterusingfdesign.lowpass.
2-53
2Filter Design and Implementation
Set the filter specifications using the 'N,Fc' specification string. With
fdesign.lowpass, you can specify your filter design in Hz.
Fs = 1000;
d = fdesign.lowpass('N,Fc',20,150,Fs);
Hd = design(d,'window','Window',kaiser(21,3));
Filter the data and plot the result.
y1 = filter(Hd,x);
plot(t(1:100),x(1:100),'k');
hold on;
plot(t(1:100),y1(1:100),'r','linewidth',2);
legend('Original Signal','Filtered Data','Location','SouthEast');
xlabel('Seconds'); ylabel('Amplitude');
Design and implement a lowpass FIR filter using the window method with
the interactive tool fdatool.
Start FDATool by entering
fdatool
at the command line.
Set the Response Type to Lowpass.SettheDesign Method to FIR and
select the Window method.
Under Filter Order,selectSpecify order. Set the order to 20.
Under Frequency Specifications.SetUnits to Hz,Fs: to 1000, and Fc: to
150.
2-54
Filtering Data With Signal Processing Toolbox™ Software
Click Design Filter.
Select File —>Export... to export your FIR filter to the MATLAB workspace
as coefficients or a filter object. In this example, export the filter as an object.
Specify the variable name as Hd1.
2-55
2Filter Design and Implementation
Click Export.
Filter the input signal in the command window with the exported filter object.
Plot the result for the first ten periods of the 100-Hz sinusoid.
y2 = filter(Hd1,x);
plot(t(1:100),x(1:100),'k');
hold on;
plot(t(1:100),y1(1:100),'r','linewidth',2);
legend('Original Signal','Filtered Data','Location','SouthEast');
xlabel('Seconds'); ylabel('Amplitude');
Select File —> Generate MATLAB Code to generate a MATLAB function
to create a filter object using your specifications.
You can also use the interactive tool filterbuilder to design your filter.
Bandpass Filters — Minimum-Order FIR and IIR
Systems
This example shows you how to design a bandpass filter and filter data with
minimum-order FIR equiripple and IIR Butterworth filters. The example uses
fdesign.bandpass and the interactive tool fdatool.
You can model many real-world signals as a superposition of oscillating
components, a low-frequency trend, and additive noise. For example,
2-56
Filtering Data With Signal Processing Toolbox™ Software
economic data often contain oscillations, which represent cycles superimposed
on a slowly varying upward or downward trend. In addition, there is an
additive noise component, which is a combination of measurement error and
the inherent random fluctuations in the process.
In these examples, assume you sample some process every day for 1 year.
Assume the process has oscillations on approximately one-week and
one-month scales. In addition, there is a low-frequency upward trend in the
data and additive N(0,1/4) white Gaussian noise.
Create the signal as a superposition of two sine waves with frequencies of 1/7
and 1/30 cycles/day. Add a low-frequency increasing trend term and N(0,1/4)
white Gaussian noise. Set the random number generator to the default state
for reproducible results. The data is sampled at 1 sample/day. Plot the
resulting signal and the power spectral density (PSD) estimate.
rng default;
Fs =1;
n = 1:365;
x = cos(2*pi*(1/7)*n)+cos(2*pi*(1/30)*n-pi/4);
trend = 3*sin(2*pi*(1/1480)*n);
y = x+trend+0.5*randn(size(n));
subplot(211);
plot(n,y); xlabel('Days'); set(gca,'xlim',[1 365]);
grid on;
subplot(212);
psdest = psd(spectrum.periodogram,y,'Fs',Fs);
plot(psdest.Frequencies,10*log10(psdest.Data));
xlabel('Cycles/day'); ylabel('dB'); grid on;
2-57
2Filter Design and Implementation
The low-frequency trend appears in the power spectral density estimate
as increased low-frequency power. The low-frequency power appears
approximately 10 dB above the oscillation at 1/30 cycles/day. Use this
information in the specifications for the filter stopbands.
Design minimum-order FIR equiripple and IIR Butterworth filters with the
following specifications: passband from [1/40,1/4] cycles/day and stopbands
from [0,1/60] and [1/4,1/2] cycles/day. Set both stopband attenuations to 10 dB
and the passband ripple tolerance to 1 dB.
d = fdesign.bandpass('Fst1,Fp1,Fp2,Fst2,Ast1,Ap,Ast2',...
1/60,1/40,1/4,1/2,10,1,10,1);
Hd1 = design(d,'equiripple');
Hd2 = design(d,'butter');
2-58
Filtering Data With Signal Processing Toolbox™ Software
Compare the order of the FIR and IIR filters and the unwrapped phase
responses.
fprintf('The order of the FIR filter is %d\n',length(Hd1.Numerator)-1);
[b,a] = sos2tf(Hd2.sosMatrix,Hd2.ScaleValues);
fprintf('The order of the IIR filter is %d\n',length(max(b,a))-1);
[phifir,w] = phasez(Hd1,[],1);
[phiiir,w] = phasez(Hd2,[],1);
plot(w,unwrap(phifir),'b'); hold on;
plot(w,unwrap(phiiir),'r'); grid on;
xlabel('Cycles/Day'); ylabel('Radians');
legend('FIR Equiripple Filter','IIR Butterworth Filter');
2-59
2Filter Design and Implementation
The IIR filter has a much lower order that the FIR filter. However, the FIR
filter has a linear phase response over the passband, while the IIR filter does
not. The FIR filter delays all frequencies in the filter passband equally, while
the IIR filter does not.
Additionally, the rate of change of the phase per unit of frequency is greater
in the FIR filter than in the IIR filter.
Design a lowpass FIR equiripple filter for comparison. The lowpass filter
specifications are: passband [0,1/4] cycles/day, stopband attenuation equal to
10 dB, and the passband ripple tolerance set to 1 dB.
dlow = fdesign.lowpass('Fp,Fst,Ap,Ast',1/4,1/2,1,10,1);
Hdlow = design(dlow,'equiripple');
Filter the data with the bandpass and lowpass filters.
yfir = filter(Hd1,y);
yiir = filter(Hd2,y);
ylow = filter(Hdlow,y);
Plot the PSD estimate of the bandpass IIR filter output. You can replace
yiir with yfir in the following code to view the PSD estimate of the FIR
bandpass filter output.
psdest = psd(spectrum.periodogram,yiir,'Fs',Fs);
plot(psdest.Frequencies,10*log10(psdest.Data));
xlabel('Cycles/day'); ylabel('dB'); grid on;
2-60
Filtering Data With Signal Processing Toolbox™ Software
The PSD estimate shows the bandpass filter attenuates the low-frequency
trend and high-frequency noise.
Plot the first 120 days of FIR and IIR filter output.
plot(n(1:120),yfir(1:120),'b');
hold on;
plot(n(1:120),yiir(1:120),'r');
xlabel('Days'); axis([1 120 -2.8 2.8]);
legend('FIR bandpass filter output','IIR bandpass filter output',...
'Location','SouthEast');
2-61
2Filter Design and Implementation
The increased phase delay in the FIR filter is evident in the filter output.
Plot the lowpass FIR filter output superimposed on the superposition of the
7-day and 30-day cycles for comparison.
plot(n,x,'k');
hold on;
plot(n,ylow,'r'); set(gca,'xlim',[1 365]);
legend('7-day and 30-day cycles','FIR lowpass filter output',...
'Location','NorthWest');
xlabel('Days');
2-62
Filtering Data With Signal Processing Toolbox™ Software
You can see in the preceding plot that the low-frequency trend is evident in
the lowpass filter output. While the lowpass filter preserves the 7-day and
30-day cycles, the bandpass filters perform better in this example because the
bandpass filters also remove the low-frequency trend.
Design and implement the bandpass Butterworth (IIR) filter with the
interactive tool fdatool.
Start FDATool by entering
fdatool
at the command line.
2-63
2Filter Design and Implementation
Set the Response Type to Bandpass.SettheDesign Method to IIR and
select the Butterworth design.
Under Filter Order, select Minimum order.
Under Frequency Specifications.SetUnits to Hz,Fs: to 1 , Fstop1: to
1/60, Fpass1: to 1/40, Fpass2: to 1/4, and Fstop2: to 1/2. Under Magnitude
Specifications,setAstop1: and Astop2: to 10 and Apass: to 1.
2-64
Filtering Data With Signal Processing Toolbox™ Software
Click Design Filter.
Select File —> Export... to export your IIR filter to the MATLAB workspace
as coefficients or a filter object. In this example, export the filter as an object.
Specify the variable name as Hd3.
Click Export.
Filter the input signal in the command window with the exported filter object.
yfilt = filter(Hd3,x);
SelectFile —> Generate MATLAB Code to generate a MATLAB function
to create a filter object using your specifications.
You can also use the interactive tool filterbuilder to design your filter.
Zero-Phase Filtering
These examples show you how to perform zero-phase filtering. The signal and
filters are described in “Lowpass FIR Filter — Window Method” on page 2-52
and “Bandpass Filters — Minimum-Order FIR and IIR Systems” on page 2-56.
2-65
2Filter Design and Implementation
Repeat the signal generation and lowpass filter design with fir1 and
fdesign.lowpass. You do not have to execute the following code if you
already have these variables in your workspace.
rng default;
Fs = 1000;
t = linspace(0,1,Fs);
x = cos(2*pi*100*t)+0.5*randn(size(t));
% Using fir1
fc = 150;
Wn = (2/Fs)*fc;
b = fir1(20,Wn,'low',kaiser(21,3));
% Using fdesign.lowpass
d = fdesign.lowpass('N,Fc',20,150,Fs);
Hd = design(d,'window','Window',kaiser(21,3));
Filter the data using filter. Plot the first 100 points of the filter output
along with a superimposed sinusoid with the same amplitude and initial
phase as the input signal.
yout = filter(Hd,x);
xin = cos(2*pi*100*t);
plot(t(1:100),xin(1:100),'k');
hold on; grid on;
plot(t(1:100),yout(1:100),'r','linewidth',2);
xlabel('Seconds'); ylabel('Amplitude');
legend('Input Sine Wave','Filtered Data',...
'Location','NorthEast');
2-66
Filtering Data With Signal Processing Toolbox™ Software
Looking at the initial 0.01 seconds of the filtered data, you see that the output
is delayed with respect to the input. The delay appears to be approximately
0.01 seconds, which is almost 1/2 thelengthoftheFIRfilterinsamples
(10*0.001).
This delay is due to the filter’s phase response. The FIR filter in these
examples is a type I linear-phase filter. The group delay of the filter is 10
samples.
Plot the group delay using fvtool.
fvtool(Hd,'analysis','grpdelay');
2-67
2Filter Design and Implementation
In many applications, phase distortion is acceptable. This is particularly true
when phase response is linear. In other applications, it is desirable to have
a filter with a zero-phase response. A zero-phase response is not technically
possibly in a noncausal filter. However, you can implement zero-phase
filtering using a causal filter with filtfilt.filtfilt does not accept filter
objects as input arguments, but you can still use filter objects with filtfilt
by extracting the filter coefficients from the object.
Filter the input signal using the coefficients in the filter object. Plot the
responses to compare the filter outputs obtained with filter and filtfilt.
yzp = filtfilt(Hd.Numerator,1,x);
% or yzp = filtfilt(b,1,x);
plot(t(1:100),xin(1:100),'k');
hold on;
plot(t(1:100),yout(1:100),'r','linewidth',2);
plot(t(1:100),yzp(1:100),'b','linewidth',2);
xlabel('Seconds'); ylabel('Amplitude');
legend('100-Hz Sine Wave','Filtered Signal','Zero-phase Filtering',...
'Location','NorthEast');
2-68
Filtering Data With Signal Processing Toolbox™ Software
In the preceding figure, you can see that the output of filtfilt does not
exhibit the delay due to the phase response of the FIR filter.
The IIR bandpass filter designed in “Bandpass Filters — Minimum-Order
FIR and IIR Systems” on page 2-56 is a biquad filter. Stated equivalently,
the IIR filter is in the form of cascaded second-order sections. To implement
zero-phase filtering with a discrete-time biquad filter, you must input the
matrix of second-order sections and the gain values for each of those sections
into filtfilt.
Zero phase filter the data in “Bandpass Filters — Minimum-Order FIR and
IIR Systems” on page 2-56 with the IIR bandpass filter. For convenience, the
2-69
2Filter Design and Implementation
code to generate the signal and filter is repeated. You do not have to execute
this code if you already have these variables in your workspace.
Generate the data.
rng default;
Fs =1;
n = 1:365;
x = cos(2*pi*(1/7)*n)+cos(2*pi*(1/30)*n-pi/4);
trend = 3*sin(2*pi*(1/1480)*n);
y = x+trend+0.5*randn(size(n));
Specify and design the filter.
d = fdesign.bandpass('Fst1,Fp1,Fp2,Fst2,Ast1,Ap,Ast2',...
1/60,1/40,1/4,1/2,10,1,10,1);
Hd2 = design(d,'butter');
Use filtfilt to zero-phase filter the input. Input the matrix of second-order
sections and gain (scale) values along with your signal.
yzpiir = filtfilt(Hd2.sosMatrix,Hd2.ScaleValues,y);
2-70
Selected Bibliography
Selected Bibliography
[1] Karam, L.J., and J.H. McClellan. “Complex Chebyshev Approximation for
FIR Filter Design.” IEEE Trans. on Circuits and Systems II. March 1995.
[2] Selesnick, I.W., and C.S. Burrus. “Generalized Digital Butterworth
Filter Design.” Proceedings of the IEEE Int. Conf. Acoust., Speech, Signal
Processing. Vol. 3 (May 1996).
[3] Selesnick, I.W., M. Lang, and C.S. Burrus. “Constrained Least Square
Design of FIR Filters without Specified Transition Bands.” Proceedings of
the IEEE Int. Conf. Acoust., Speech, Signal Processing. Vol. 2 (May 1995).
Pgs. 1260-1263.
2-71
2Filter Design and Implementation
2-72
3
Designing a Filter in
Fdesign — Process
Overview
3Designing a Filter in Fdesign — Process Overview
Process Flow Diagram and Filter Design Methodology
In this section...
“Exploring the Process Flow Diagram” on page 3-2
“Selecting a Response” on page 3-4
“Selecting a Specification” on page 3-4
“Selecting an Algorithm” on page 3-6
“Customizing the Algorithm” on page 3-8
“Designing the Filter” on page 3-8
“Design Analysis” on page 3-9
“Realize or Apply the Filter to Input Data” on page 3-10
Note You must minimally have the Signal Processing Toolbox installed
to use fdesign and design. Some of the features described below may
be unavailable if your installation does not additionally include the DSP
System Toolbox™ license. The DSP System Toolbox significantly expands the
functionality available for the specification, design, and analysis of filters.
You can verify the presence of both toolboxes by typing ver at the command
prompt.
Exploring the Process Flow Diagram
The process flow diagram shown in the following figure lists the steps and
shows the order of the filter design process.
3-2
Process Flow Diagram and Filter Design Methodology
The first four steps of the filter design process relate to the filter Specifications
Object, while the last two steps involve the filter Implementation Object. Both
of these objects are discussed in more detail in the following sections. Step 5
- the design of the filter, is the transition step from the filter Specifications
Object to the Implementation object. The analysis and verification step is
3-3
3Designing a Filter in Fdesign — Process Overview
completely optional. It provides methods for the filter designer to ensure that
the filter complies with all design criteria. Depending on the results of this
verification, you can loop back to steps 3 and 4, to either choose a different
algorithm, or to customize the current one. You may also wish to go back to
steps 3 or 4 after you filter the input data with the designed filter (step 7),
andfindthatyouwishtotweakthe filter or change it further.
The diagram shows the help command for each step. Enter the help line at the
MATLAB command prompt to receive instructions and further documentation
links for the particular step. Not all of the steps have to be executed explicitly.
For example, you could go from step 1 directly to step 5, and the interim three
steps are done for you by the software.
The following are the details for each of the steps shown above.
Selecting a Response
If you type:
help fdesign/responses
at the MATLAB command prompt, you see a list of all available filter
responses. The responses marked with an asterisk require the DSP System
Toolbox.
You must select a response to initiate the filter. In this example, a bandpass
filter Specifications Object is created by typing the following:
d = fdesign.bandpass
Selecting a Specification
Aspecification is an array of design parameters for a given filter. The
specification is a property of the Specifications Object.
Note A specification is not the same as the Specifications Object. A
Specifications Object contains a specification as one of its properties.
3-4
Process Flow Diagram and Filter Design Methodology
When you select a filter response, there are a number of different
specifications available. Each one contains a different combination of design
parameters. After you create a filter Specifications Object, you can query the
available specifications for that response. Specifications marked with an
asterisk require the DSP System Toolbox.
>> d = fdesign.bandpass; % step 1 - choose the response
>> set (d, 'specification')
ans =
'Fst1,Fp1,Fp2,Fst2,Ast1,Ap,Ast2'
'N,F3dB1,F3dB2'
'N,F3dB1,F3dB2,Ap'
'N,F3dB1,F3dB2,Ast'
'N,F3dB1,F3dB2,Ast1,Ap,Ast2'
'N,F3dB1,F3dB2,BWp'
'N,F3dB1,F3dB2,BWst'
'N,Fc1,Fc2'
'N,Fp1,Fp2,Ap'
'N,Fp1,Fp2,Ast1,Ap,Ast2'
'N,Fst1,Fp1,Fp2,Fst2'
'N,Fst1,Fp1,Fp2,Fst2,Ap'
'N,Fst1,Fst2,Ast'
'Nb,Na,Fst1,Fp1,Fp2,Fst2'
>> d=fdesign.arbmag;
>> set(d,'specification')
ans =
'N,F,A'
'N,B,F,A'
The set command can be used to select one of theavailablespecificationsas
follows:
>> d = fdesign.lowpass; % step 1
>> % step 2: get a list of available specifications
>> set (d, 'specification')
3-5