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 .
Page Count: 527 [warning: Documents this large are best viewed by clicking the View PDF Link!]

toc
tables

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

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

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

Filter Design and Implementation

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

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

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

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

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

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

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

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

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

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

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

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

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

Prediction Polynomial ............................. 13-2

Formant Estimation with LPC Coefficients .......... 13-5

AR Order Selection with Partial Autocorrelation

Sequence ....................................... 13-9

Transforms

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

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

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

Technical Conventions

Index

xvi Contents

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

() ()( )=−

=−∞

∞

∑

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

() () () () () ... ( )

() () ...

++++

+++

−−

−

12 1

1(() ()

+−

where H(z)isthefilter’stransfer 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

() ()() ( )

() ()() ( ) ()()

=+−

=+−−

−

212



221

111() ( )() ( ) ( )()

() ()() ()

mbn xmz m an ym

zmbnxman

=− + −−−

=−

−

−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

≠

⎧

⎨

⎩

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

() ( ) ( ) ... ( )

() () ... (

=++++

++++

−−

−

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

 



() ()

=−

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

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

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

() () () ( )

() () ( ) ()=++…++

++…++

−−

12 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

() ()

()

( ( ))( ( ))...( ( ))

( ( ))( ( ))

−− −

−−

....( ( ))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

()=++

+++

=++

+++

−−

−−−

23 4

13 3

234

331

123

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

( ) () ()

() () ()

+= +

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

pnz

kkz

()

() ... ()

() ( ) ( ) ...=−++

−++ ++

−−

−

11 1 12

1kkm n z mn

()

−+ −−

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:

pjz

rj s

pjz

()

(())

... ()

(())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

()=−+

−

−−

16 8

1-26

Linear System Models

b = [-4 8];

a=[168];

[r,p,k] = residuez(b,a)

which corresponds to

()=−

+++

−−

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

() ( ) ( ) ,,,=− − + − = …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

Lkk k

kk k

() ()==

−−

∏∏

122

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

01 11 21 01 11 21

02 12 22 02 12 22

......

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

bs bs bn

as as a

() ()

()

() () ( )

() () (

== ++…++

++…+

−

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

() ()

()

( ( ))( ( )) ( ( ))

( ( ))( ( )) (

==−−…−

−−…

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

Nkn

() (),+= +

−

∑

and

xn NXk W

Nkn

() .() 







111

In these equations, the series subscripts begin with 1 instead of 0 because of

the MATLAB vector indexing scheme, and

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

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)

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)

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)

0.0317 0.0951 0.0951 0.0317

1.0000 -1.4590 0.9104 -0.1978

[b,a] = butter(3,0.25)

0.0317 0.0951 0.0951 0.0317

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)

0.0950 0.2849 0.2849 0.0950

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

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

restriction

Type II Odd even:

bk bn k k n( ) ( ), ,...,=+− = +211

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

jn jn

() ( ) sin( )

===

−−

∫∫

πωω

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

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/

[numt,dent] = lp2lp (num,den,Wo)

[At,Bt,Ct,Dt] =lp2lp (A,B,C,D,Wo)

Lowpass to highpass

′=ss

[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

′=+

ωω

1(/ )

[numt,dent] = lp2bp (num,den,Wo,Bw)

[At,Bt,Ct,Dt] =lp2bp (A,B,C,D,Wo,Bw)

Lowpass to bandstop

′=+

sBs

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

ωω

=−

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

() ()==−

Bilinear transformation maps the jΩ-axis of the continuous domain to the

unit circle of the discrete domain according to

=⎛

⎝

⎜⎞

⎠

⎟

−

tan Ω

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

[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

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

3Designing a Filter in Fdesign — Process Overview

ans =

'Fp,Fst,Ap,Ast'

'N,F3dB'

'N,F3dB,Ap'

'N,F3dB,Ap,Ast'

'N,F3dB,Ast'

'N,F3dB,Fst'

'N,Fc'

'N,Fc,Ap,Ast'

'N,Fp,Ap'

'N,Fp,Ap,Ast'

'N,Fp,F3dB'

'N,Fp,Fst'

'N,Fp,Fst,Ap'

'N,Fp,Fst,Ast'

'N,Fst,Ap,Ast'

'N,Fst,Ast'

'Nb,Na,Fp,Fst'

>> %step 2: set the required specification

>> set (d, 'specification', 'N,Fc')

If you do not perform this step explicitly, fdesign returns the default

specification for the response you chose in “Select a Response” on page 4-3, and

provides default values for all design parameters included in the specification.

Selecting an Algorithm

The availability of algorithms depends the chosen filter response, the design

parameters, and the availability of the DSP System Toolbox. In other words,

for the same lowpass filter, changing the specification string also changes the

available algorithms. In the following example, for a lowpass filter and a

specification of 'N, Fc', only one algorithm is available—window.

>> %step 2: set the required specification

>> set (d, 'specification', 'N,Fc')

>> designmethods (d) %step3: get available algorithms

3-6

Process Flow Diagram and Filter Design Methodology

Design Methods for class fdesign.lowpass (N,Fc):

window

However, for a specification of 'Fp,Fst,Ap,Ast',anumberofalgorithmsare

available. If the user has only the Signal Processing Toolbox installed, the

following algorithms are available:

>>set (d, 'specification', 'Fp,Fst,Ap,Ast')

>>designmethods(d)

Design Methods for class fdesign.lowpass (Fp,Fst,Ap,Ast):

butter

cheby1

cheby2

ellip

equiripple

kaiserwin

If the user additionally has the DSP System Toolbox installed, the number of

available algorithms for this response and specification string increases:

>>set(d,'specification','Fp,Fst,Ap,Ast')

>>designmethods(d)

Design Methods for class fdesign.lowpass (Fp,Fst,Ap,Ast):

butter

cheby1

cheby2

ellip

equiripple

ifir

kaiserwin

3-7