Mathworks SIGNAL PROCESSING TOOLBOX 6 user guide

Download

Signal Processing

User’s Guide

Toolbox™ 6

How to Contact The MathWorks

www.mathworks. comp.soft-sys.matlab Newsgroup www.mathworks.com/contact_TS.html T echnical Support

suggest@mathworks.com Product enhancement suggestions bugs@mathwo doc@mathworks.com Documentation error reports service@mathworks.com Order status, license renewals, passcodes info@mathwo

com

rks.com

Web

Bug reports

Sales, prici

ng, 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 MathW orks, 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 Docu mentation, 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

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

Contents

Filtering, Linear Systems and Transforms

Overview

Filter Implem

Filtering Ove Convolution Filters and T Filtering wi

The filter F

Other Func

Multirate Anti-Caus Frequenc

Impulse

Frequen

Digita Analog Magnit Delay

Pole Analysis

Zero-

unction

tions for Filtering

Filter Bank Implementation al, Zero-Phase Filter Implementation

y Domain Filter Implementation

Response

cy Response

lDomain

Domain

ude and Phase

...........................................

entation and Analysis

rview

and Filtering

ransfer Functions

th the filter Function

................................

...........................

.................

......................

.....................

................................

.......................

................

.............

.................................

...............................

...................................

..............................

.................................

.........

1-2 1-2 1-2 1-3 1-4

1-6

1-8 1-8 1-9

1-10

1-12

1-14 1-14 1-16 1-17 1-18

1-21

ar System Models

Line

lable Models

Avai

rete-Time System Models

Disc

tinuous-Time System Models

Con

ear System Transformations

Lin

screte Fourier Transform

..................................

..............................

.......................

.....................

........................

1-23 1-23 1-23 1-3 1-3

2 3

Filter Design and Implementation

Filter Requirements and Specification .............. 2-2

IIR Filter Design

IIR vs. FIR Filters Classical IIR Filters Other IIR Filters IIR Filter Method Summary Classical IIR Filter Design Using Analog Prototyping Comparison of Classical IIR Filter Types

FIR Filter Design

FIR vs. IIR Filters FIR Filter Summary Linear Phase Filters Windowing Method Multiband FIR Filter Design with Transition Bands Constrained Least Squares FIR Filter Design Arbitrary-Response Filter Design

Special Topics in IIR Filter Design

Classic IIR Filter Design Analog Proto typ e Design Frequency Transformation Filter Discretization

.................................. 2-4

................................. 2-4

............................... 2-4

.................................. 2-5

......................... 2-5

.............. 2-9

.................................. 2-17

................................. 2-17

............................... 2-18

................................ 2-20

.................... 2-37

.................. 2-43

........................... 2-43

........................... 2-44

.......................... 2-44

............................... 2-46

.... 2-6

..... 2-24

.......... 2-31

vi Contents

Selected Bibliography

.............................. 2-52

Designing a Filter in Fdesign — Process

Process Flow Diagram and Filter Design

Methodology

Exploring the Process Flow Diagram Selecting a Response

.................................... 3-2

............................... 3-4

Overview

.................. 3-2

Selecting a Specification ............................ 3-4

Selecting an Algorithm Customizing the Algorithm Designing the Filter Design Analysis RealizeorApplytheFiltertoInputData

................................... 3-9

............................. 3-6

......................... 3-8

............................... 3-8

.............. 3-10

Designing a Filter in the Filterbuilder GUI

The Graphical Interface to Fdesign ................. 4-2

Introduction to Filterbuilder Filterbuilder Design Process Select a Response Select a Specification Select an Algorithm Customize the Algorithm Analyze the Design RealizeorApplytheFiltertoInputData

................................. 4-3

.............................. 4-5

................................ 4-5

................................ 4-8

........................ 4-2

........................... 4-6

.............. 4-8

Designing a FIR Filter Using filterbuilder

FIR Filter Design

................................. 4-10

........... 4-10

FDATool: A Filter Design and Analysis GUI

Overview ......................................... 5-2

Introduction to FDA Tool Integrated Products Filter Design Methods Using the Filter Design and Analysis Tool Analyzing Filter Responses Filter Design and Analysis Tool Panels Getting Help

Opening FDATool

..................................... 5-6

.................................. 5-7

........................... 5-2

............................... 5-2

............................. 5-3

............. 5-4

......................... 5-5

................ 5-5

vii

Choosing a Response Type ......................... 5-8

Choosing a Filter Design Method

Setting the Filter Design Specifications

Viewing Filter Specifications Filter Order Options Bandpass Filter Frequency Specifications Bandpass Filter Magnitude Specifications

Computing the Filter Coefficients

Analyzing the Filter

Displaying Filte r Responses Using Data Tips Drawing Spectral Masks Changing the Sampling Frequency DisplayingtheResponseinFVTool

Editing the Filter Using the Pole/Zero Editor

Displaying the Pole-Zero Plot Changing the Pole-Zero Plot

...................................... 5-10

.......................................... 5-11

............................... 5-16

.................................. 5-18

............................ 5-19

.................... 5-9

............. 5-10

........................ 5-10

.............. 5-12

............. 5-13

................... 5-15

......................... 5-16

................... 5-20

................... 5-21

........................ 5-23

......................... 5-24

........ 5-23

viii Contents

Converting the Filter Structure

Converting to a New Structure Converting to Second-Order Sections

Importing a Filter Design

Import Filter Panel Filter Structures

Exporting a Filter Design

Exporting Coefficients or Objects to the Workspace Exporting Coefficients to an ASCII File Exporting Coefficients or Objects to a MAT-File Exporting to SPTool Exporting to a Simulink Model Other Ways to Export a Filter

................................ 5-30

.................................. 5-32

............................... 5-37

..................... 5-27

...................... 5-27

................. 5-28

.......................... 5-30

.......................... 5-35

............... 5-36

...................... 5-38

....................... 5-41

...... 5-35

......... 5-37

Generating a C Header File ......................... 5-42

Generating an M-File

Managing Filters in the Current Session

Saving and Opening Filter Design Sessions

.............................. 5-44

............. 5-45

.......... 5-47

Statistical Signal Processing

Correlation and Covariance ........................ 6-2

Background Information Using xcorr and xcov Functions Bias and Normalization Multiple Channels

Spectral Analysis

Background Information Spectral Estimation Method Nonparametric Methods Parametric Methods Using FFT to Obtain Simple Spectral Analysis Plots

.................................. 6-5

............................ 6-2

...................... 6-3

............................ 6-3

................................. 6-4

............................ 6-5

......................... 6-7

............................ 6-9

............................... 6-32

..... 6-45

Selected Bibliography

.............................. 6-48

Special Topics

Windows .......................................... 7-2

Why Use Windows? Available Window Functions Graphical User Interface Tools Basic Shapes

................................ 7-2

........................ 7-2

...................... 7-3

..................................... 7-3

Generalized Cosine Windows ........................ 7-7

Kaiser Window Chebyshev Window

................................... 7-9

................................ 7-14

Parametric Modeling

What is Parametric Modeling Available Parametric Modeling Functions Time-Domain Based Modeling Frequency-Domain Based Modeling

Resampling

Available Resampling Functions resample Function decimate and interp Functions upfirdn Function spline Function

Cepstrum Analysis

What Is a Cepstrum? Inverse Com plex Cepstrum

FFT-Based Time-Frequency Analysis

Median Filtering

Communications A pplications

Modulation Demodulation Voltage Controlled Oscillator

....................................... 7-25

....................................... 7-34

.................................... 7-35

.............................. 7-15

........................ 7-15

....................... 7-16

..................... 7-25

................................. 7-25

....................... 7-27

.................................. 7-27

................................... 7-27

................................. 7-28

.............................. 7-28

......................... 7-31

.................................. 7-33

...................... 7-34

........................ 7-38

............. 7-15

.................. 7-21

................ 7-32

x Contents

Deconvolution

Specialized Transforms

Chirp z-Trans form Discrete Cosine Transform Hilbert Transform Walsh–Hadamard Transform

Selected Bibliography

..................................... 7-39

................................. 7-40

................................. 7-44

............................ 7-40

.......................... 7-41

........................ 7-45

.............................. 7-51

SPTool: A Signal Processing GUI Suite

SPTool: An Interactive Signal Processing

Environment

SPTool Ove rview SPTool Da ta Structures

.................................... 8-2

.................................. 8-2

............................ 8-3

Opening SPTool

Getting Context-Sensitive Help

Signal Browser

Overview of the Signal Browser Opening the Signal Browser

FDATool

Filter Visualization Tool

Connection between FVTool and SPTool Opening the Filter Visualization Tool Analysis Parameters

Spectrum Viewer

Spectrum Viewer Overview Opening the Spectrum Viewer

Filtering and Analysis of Noise

Overview Importing a Signal into SPTool Designing a Filter Applying a Filter to a Signal Analyzing a Signal Spectral Analysis in the Spectrum Viewer

.......................................... 8-10

................................... 8-4

..................... 8-6

.................................... 8-7

...................... 8-7

......................... 8-7

........................... 8-12

............................... 8-13

.................................. 8-14

......................... 8-14

....................... 8-14

..................... 8-17

........................................ 8-17

...................... 8-17

................................. 8-19

........................ 8-21

................................ 8-23

............... 8-12

................. 8-12

............. 8-25

Exporting Signals, Filters, and Spectra

Opening the Export Dialog Box Exporting a Filter to the MATLAB Workspace

...................... 8-28

.............. 8-28

.......... 8-29

Accessing Filter Parameters ........................ 8-30

Accessing Filter Parameters in a Saved Filter Accessing Parameters in a Saved Spectrum

.......... 8-30

............ 8-31

Importing Filters and Spectra

Similarities to Other Procedures Importing Filters Importing Spectra

Loading Variables from the Disk

Saving and Loading Sessions

SPTool Sessions Filter Formats

Selecting Signals, Filters, and Spectra

Editing Signals, Filters, or Spectra

Making Signal Measurements with Markers

Setting Preferences

Overview o f Setting Preferences Summary of Settable Preferences Setting the Filter Design Tool

.................................. 8-33

................................. 8-35

................................... 8-38

.................................... 8-38

................................ 8-44

...................... 8-33

..................... 8-33

.................... 8-37

....................... 8-38

............... 8-40

.................. 8-41

..................... 8-44

.................... 8-45

....................... 8-46

......... 8-42

xii Contents

Using the Filter Designer

Filter Designer Filter Types FIR Filter Methods IIR Filter Methods Pole/Zero Editor Spectral Overlay Feature Opening the Filter Designer Accessing Filter Parameters in a Saved Filter Designing a Filter with the Pole/Zero Editor Positioning Poles and Zeros Redesigning a Filter Using the Magnitude Pl ot

................................... 8-48

...................................... 8-48

................................... 8-49

........................... 8-48

................................ 8-49

........................... 8-49

......................... 8-49

......................... 8-55

.......... 8-51

........... 8-54

......... 8-57

Embedded MATLAB Support in Signal

Processing Toolbox

Supported Functions .............................. 9-2

Specifying Inputs in Embedded MATLAB

Defining Input Size and Type Inputs must be Constants

Embedded M ATL AB Examples

Apply Window to Input Signal Apply Lowpass Filter to Input Signal Cross Correlate or Autocorrelate Input Data

freqz W ith No Output Arguments ................... 9-25

Zero Phase Filtering

............................... 9-26

........................ 9-16

........................... 9-17

...................... 9-20

....................... 9-20

................. 9-23

........... 9-16

........... 9-23

Function Reference

Signal Processing Functions in MATLAB ............ 10-3

Digital Filters

Discrete-Time Filters FIR Filter Design Communications Filters IIR Digital Filter Design IIR Filter Order Estimation Filter Analysis Filter Implementation Filter Specification Objects -- Response Types Filter Specification Objects -- Design Methods

..................................... 10-3

.............................. 10-4

................................. 10-5

............................ 10-6

......................... 10-6

.................................... 10-7

............................. 10-7

.......... 10-8

.......... 10-9

Analog Filters

Analog Low p ass Filter Prototypes Analog Filter Design Filter Analysis Analog Filter Transformation

..................................... 10-9

............................... 10-10

.................................... 10-10

.................... 10-9

....................... 10-10

xiii

Filter Discretization ............................... 10-11

Linear Systems

Windows

Transforms

Cepstral Analysis

Statistical Signal Processing

Parametric Modeling

Linear Prediction

Multirate Signal Processing

Waveform Generation

Specialized Operations

.......................................... 10-12

.................................... 10-11

....................................... 10-13

.................................. 10-14

.............................. 10-16

.................................. 10-16

.............................. 10-18

............................. 10-19

........................ 10-14

........................ 10-18

xiv Contents

GUIs

.............................................. 10-20

By Category

Windows .......................................... 11-2

Functions — Alphabetical List

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-35

1 Filtering, 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 F unction” on page 1-4

Filtering Overview

This section d escribes 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.

1-2

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 tw o- dimensional signal processing using the

A digital filter’s output y(k)isrelatedtoitsinputx(k) by convolution with its impulse response h ( k).

yk hlxk l

conv2 function.

∞

() ()( )=−

∑

=−∞

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 vector

x, h(k)inavectorh, a nd 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);

conv.Storex(k)ina

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-transfo rm X(z)oftheinputby

()

m+−

Yz HzXz

() () ()

−−

12 1

bbz bnz

() () ... ( )

++++

aaz a

() () ...

+++

−

(()

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 n and 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 a and b coefficients present:

• When

n =0(thatis,b is a scalar), the filter is an Infinite Impulse Re sp on s e

(IIR), all-pole, recursive, or autoregressive (AR) filter.

m =0(thatis,a is a scalar), the filter is a Finite Impulse Response

(FIR), all-zero, nonrecursive, or moving-average (MA) filter.

1-3

1 Filtering, Linear Systems and Transforms Overview

• If both n and m are 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 sto chastic 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 pro gressio n is

1-4

ybx

() () ()

111

ybxbxay

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

2122121

=+−

ybx

() ()()

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 B and A represent the coefficients of a filter in transfer function form. Note that the and input terms are separated in the difference equation. F o r the example, the previous coefficient vectors represent a linear constant-coefficient difference equation of

A coefficient vectors are written a s if the output

Filter Implementation and Analysis

yn yn xn() . ( ) ()−−=09 1

Changing the sign of the A(2) coefficie nt, 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

is, the length of

gives you as many output samples as there are input samples, that

y is the same as the length o f x.Ifthefirstelementofa

is not 1, filter divides the coefficients by a(1) before implementing the difference equation.

1-5

1 Filtering, 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

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

=+−−



−

zmbnxman

−

212

mbn xmz m an ym

=− + −−−

221

() ()() ()

=−

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

−

yym()

1-6

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

The filter Function

x = [x1;x2]. If there is not sufficient memory to hold the combined sequence,

filter the subsequences the filtered sequences, use the final conditions from to filter

x2:

[y1,zf] = filter(b,a,x1); y2 = filter(b,a,x2,zf);

x1 and x2 one at a time. To ensure continuity of

x1 as initial conditions

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 using

filtic,use

zf = filtic(b,a,flipud(y1),flipud(x1));

zf as above

This can be useful when filtering short data sequences, as appropriate initial conditions help reduce transient startup effects.

1-7

1 Filtering, 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 D omain 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 cascad e of three systems that performs the following tasks:

• Upsampling (zero insertion) by integer factor

• 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

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 o f 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

p/q.Set

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 th rou g h each FIR filter, resulting in a matrix of output signals.

Other functions that perform multirate filtering (with fixed fil ter) include

resample, interp,anddecimate.

Anti-Causal, Zero-Phase Filter Implementation

In the ca s e of FIR filters, it is possible to design linear phase filters that, when applied to data (using fixed number of samples. For IIR filters, however, the phase distortion is usually highly nonlinear. The signal at points before and after the current point, in essence “looking into the future,” to eliminate phase distortion.

filter or conv), simply delay the output by a

filtfilt function uses the information in the

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 x that 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 x using both filter and filtfilt for compari son :

1-9

1 Filtering, 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,':')

1-10

Both filtered versions eliminate the 40 Hz sinusoid evident in the original, solid line. The plot also shows how (

filtfilt) line is in phase with the original 3 Hz sinusoid, while the dotted

(

filter) line is delayed by about five sampl es. Also, the amplitude of the

dashed line is smaller due to the magnitude squared effects of

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.

filter and filtfilt differ; the dashed

filtfilt.

Frequency Domain Filter Implementation

Dualitybetweenthetimedomainandthefrequencydomainmakesitpossible to perform any operation in either domain. U sually one domain or the other is more convenient for a particular operation, but you can always accomplish a given operation in either domain.

Other Functions for Filtering

To implement general IIR filtering in the frequency domain, multiply the discrete Fourier transform (DFT) of the input sequ enc e 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

n increases.

For FIR filters, however, it is possible t o 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

1 Filtering, 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

⎧



()n

⎨

≠

⎩

You can gener way is

imp = [1; zeros(49,1)];

The impuls

h = filter(b,a,imp);

Asimplew Tool (

fvt

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:

ate an impulse sequence a number of ways; one straightforward

e response of the simple filter

ay to display the impulse response is with the Filter Visualization

ool

b = 1 and a = [1 -0.9] is

1-12

Impulse Response

1-13

1 Filtering, 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. Sp eci fi cal ly , the statement

[h,w] = freqz(b,a,p)

returns the p-point complex frequency response, H(ejω), of the digital fi lter.

1-14

jw jwn

()

12 1

bbe bne

( ) ( ) ... ( )

aae am

() () ... (

++++

−−

++++

−

11)e

jwm−

In its simplest form, freqz accepts the filter coefficient vectors b and a,andan integer response. actual frequency points in vector

freqz can accept other parameters, such as a sampling frequency or a

p specifying the number of points at which to calculate the frequency

freqz returns the complex frequency response in vector h,andthe

w in rad/s.

vector of arbitrary frequency points. The example below finds the 256-point frequency response for a 12th-order C hebyshev Type I filter. The call to

freqz

specifies a sampling fr equency fs of 1000 Hz:

a] = cheby1(12,0.5,200/500);

[b,

f] = freqz(b,a,256,1000);

[h,

Because the parameter list includes a sampling frequency, freqz returns a vector

f that contains the 256 frequency points between 0 and fs/2 used in

the frequency response calculation.

Frequency Response

Note This toolbox uses the convention that unit frequency is the Nyquist frequency, defined a s 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

1 Filtering, Linear Systems and Transforms Overview

1-16

freqz can also accept a vector of arbitrary frequency points for use in the

frequency response calculation. For example,

w = linsp h = freqz

calculates the complex frequency response at the frequency points in w for the filter defined by vectors 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.

ace(0,pi); (b,a,w);

b and a. The frequency points can range from 0 to

Analog Domain

freqs evaluates frequency response for an analog filter defined by two input

coefficient vectors, 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.

b and a. Its operation is similar to that of freqz;youcan

Frequency Response

Magnitude and Ph

MATLAB function frequency respo response; and phase of a Bu

[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.

angl

s are available to extract magnitude and phase from a

nse vector

returns the phase angle in radians. To extract the magnitude

tterworth filter:

ase

h. The function abs returns the magnitude of the

The unwrap function is also useful in frequency analysis. unwrap unwraps the p hase to m ake it continuous across 360º phase discontinuities by adding multiples of ±360°, as needed. To see how 25th-order lowpass FIR filter:

h = fir1(25,0.4);

Obtain the filter’s frequency response with freqz, and plot the phase in degrees:

unwrap is useful, design a

1-17

1 Filtering, Linear Systems and Transforms Overview

[H,f] = freqz(h,1,512,2); plot(f,angle(H)*180/pi); grid

1-18

It is difficult to distinguish the 360° jumps (an artifact of the arctangent function inside response.

unwrap eliminates the 360° jum ps:

plot(f,unwrap(angle(H))*180/pi);

or you can use phasez to see the unwrapped phase.

angle) from the 180° jumps that signify zeros in the frequency

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(e then the group delay is

()

=−





()



jω

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

b and a, evaluated at the frequencies in vector w.

The phase delay of a filter is the negative of phase divided by frequency:



=−

()





()

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

1 Filtering, 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)

0.9e

zer = -0.5; pol = 0.9*exp(j*2*pi*[-0.3 0.3]');

–j2π(0.3)

and

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

[b,a] = zp2tf(zer,pol,1); fvtool(b,a)

fvtool,

and click the Pole/Zero Plot toolbar button on the toolbar or select Analysis > Pole/Zero Plot to see the plot.

1-21

1 Filtering, Linear Systems and Transforms Overview

1-22

For a system in zero-pole form, supply column vector arguments z and p to

zplane:

zplane(z,p)

For a system in transfer function form, supply row vectors b and a as arguments to

zplane(b,a)

In this case zplane finds the roots of b and a using 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.

zplane:

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-32

“Linear System Transformations” on page 1-33

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, conve rt freely to and from most other models. This section prov ides a brief overview of supported linear system models and describes how to work with these models in the MATLAB technical computing environment.

Linear System Models

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-28

• “Lattice Structure” on page 1-29

• “Convolution Matrix” on p age 1-31

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

1 Filtering, 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

()=

()

−−

bbz bnz

() () ( )

++…++

12 1

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 n and 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 trans fer function is

()

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

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 k is a MATLAB scalar.

The

poly and roots functions convert between polynomial and zero-pole-gain

representations. For example, a simple IIR filter is

zq zq zqn

−− −

zp zp

( ( ))( ( ))

−−

....( ( ))zpn−

1-24

b=[234]; a = [1 3 3 1];

The zeros and poles of this filter are

q = roots(b) q=

-0.7500 + 1.1990i

-0.7500 - 1.1990i p = roots(a) p=

-1.0000

-1.0000 + 0.0000i

-1.0000 - 0.0000i k = b(1)/a(1) k=

Returning to the or iginal polynomials,

bb = k*poly(q) bb =

2.0000 3.0000 4.0000 aa = poly(p) aa =

1.0000 3.0000 3.0000 1.0000

Note that b and a in this case represent the transfer function:

−−

()=

23 4

13 3

+++

−−−

123

zzz

234

zzz

331

zzz

+++

Linear System Models

For b = [2 3 4],theroots function misses the zero for z equal to 0. In fact, it misses poles and zeros for z equal 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 repre sent 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 u is the input, x is the state vector, and y is the output. For single-channel systems,

A is an m-by-m matrix where m is the order of the filter,

1-25

1 Filtering, Linear Systems and Transforms Overview

B is a column vector, C is a row vector, and D is a scalar. State-space notation

is especially convenient for multichannel systems where input become vectors, and B, C,andD become matrices.

u and output y

State-space representation extends easily to the MATLAB environment.

C,andD are rectangular arrays; MATLAB functions treat them as individual

A, B,

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

Don’t be concerned if you are not familiar with the state-spa ce representation of linear systems. Some of the filter design algorithms use state-space form internally but do not require any knowledge of s tate-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

()

−

11 1

...

−−

−

()

pnz

()

−

kkz

( ) ( ) ...=

++ ++

kkm n z

()

−+

()

−−

1-26

provided H(z) has no repeated poles. Here, n is the degree of the denominator polynomial of the rational transfer function b(z)/a(z). If r is a pole of multiplicity s

pjz

−

,thenH(z)hastermsoftheform:

rj s

()

112 1

−− −

()

pjz

(())

−

()

...

(())1

−

+−

pjz

The Signal Processing Toolbox residuez function in converts transfer functions to and from the partial fraction e xpansion form. The “ of

residuez stands for z-domain, or discrete domain. residuez returns the

z”ontheend

Linear System Models

poles in a column vector p, the residues corresponding to the poles in a column vector vector

r, and any improper part of the original transfer function in a row 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

−

−+

()=

−−

16 8

b=[-48]; a=[168]; [r,p,k] = residuez(b,a) r=

-12

-4

-2 k=

[]

which corresponds to

()=

− +

−−

14812

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

1-27

1 Filtering, Linear Systems and Transforms Overview

imp = [1 0 0 0 0]; resptf = filter(b,a,imp) resptf =

-4 32 -160 704 -2944

respres = filter(r(1),[1 -p(1)],imp)+...

filter(r(2),[1 -p(2)],imp)

respres =

-4 32 -160 704 -2944

Second-Order Sections (SOS)

Any transfer function H(z) has a second-order sections representation

Hz H z

() ()==

∏∏

aazaz

bbzbz

0112

kk k

0112

kk k

−−

where L is the number of second-order sections that describe the system. The MATLAB environment represents the second-order se ction 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 se ction.

bbbaaa

sos

⎛

01 11 21 01 11 21

⎜

bbbaaa

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-33, perform pole-zero pairing, section scaling, and section ordering.

1-28

Linear System Models

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 N corresponding 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.

d IIR Lattice Filter structure diagrams

FIR an

general pole-zero IIR filter described by polynomial coefficients a

For a

, there are both lattice coefficients k(n) for the denominator a and

and b

er coefficients v(n) for the numerator b. The lattice/ladder filter may be

ladd

lemented as

imp

1-29

1 Filtering, Linear Systems and Transforms Overview

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];

This filter’s lattice (reflection coefficient) representation is

k = tf2latc(b) k=

0.3090

0.9801

0.0031

0.0081

-0.0082

1-30

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

b from

above and denominator polynomial:

a = [1 1/2 1/3];

The filter’s lattice representation is

[k,v] = tf2latc(b,a) k=

0.3750

0.3333

Linear System Models

0 0 0

0.6252

0.1212

0.9879

-0.0009

0.0072

-0.0082

Because abs(k) < 1 for 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) b=

1.0000 0.6149 0.9899 -0.0000 0.0031 -0.0082

The lattice or lattice/ladder coefficients can be used to implement the filter using the function

latcfilt.

Convolution Matrix

In signal process ing, 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 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 for a digital filter, is

convmtx generates a matrix whose

b, representing the numerator coefficients

b = [1 2 3]; x = randn(3,1); C = convmtx(b',3)

1-31

1 Filtering, Linear Systems and Transforms Overview

100 210 321 032 003

Two equivalent ways to convolve b with x are 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

1-32

• 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 u is a vector of nu inputs, x is an nx-element state vector, and y is a vector of ny outputs. In the MATLAB environment, separate rectangular arrays.

An equivalent representation of the state-space system is the Laplace transform transfer function description

A, B, C,andD are stored in

Ys HsUs() () ()=

where

Linear System Models

Hs CsI A B D() ( )=− +

−1

For single-input, single-output systems, this form is given by

()

bsasbs bs bn

()()() () ( )

++…++

12 1

as as a

() () (

++…+

−

mm +1)

Given the coefficients of a Laplace transform transfer function, residue determines the partial fraction expansion of the system. See the description of

residue in the MATLAB documentation for details.

The factored zero-pole-gain form is

()

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

sz sz szn

−−…−

sp sp

( ( ))( ( )) (

−−…

sspm− ())

Asinthediscrete-timecase,theMATLABenvironmentstorespolynomial coefficients in row vectors in descend ing 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; see Chapter 12, “Functions — Alphabetical List” for a complete description of each. 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.

1-33

1 Filtering, Linear Systems and Transforms Overview

Transfer Function

State-Space

Zero-PoleGain

Partial Fraction

Lattice Filter

SOS

ZeroTransfer Function

ss2tf ss2zp none none ss2sos none

zp2tf poly

residuez none none none none none

latc2tf none none none none none

sos2tf sos2ss sos2zp none none none

StateSpace

tf2ss tf2zp

zp2ss none none zp2sos none

Pole-

Gain

roots

Partial Fraction

residuez tf2latc none convmtx

Lattice Filter

SecondOrder Sections

Convolution Matrix

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. See the descriptions of the individual filter design functions in Chapter 12, “Functions — Alphabetical List” for more details.

1-34

Note AllSignalProcessingToolboxsecond-order section transformations apply only to digital filters.

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 analys is, a nd som e filter design and implementation functions) incorporate the FFT.

Discrete Fourier Transform

The MATLAB environment provides the functions

fft and ifft to compute

the discrete Fourier transform and its inverse, respectively. F or the input sequence x and 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

() () ,+= +

∑

and

−

() .()+= +

Xk W

∑

−

In these equations, the series subscripts begin with 1 instead of 0 because of the MATLAB vector indexing scheme, and

− 2/

Note TheMATLABconventionistouseanegativej for the fft function. This is an engineering convention; physics and pure mathematics typically use a positive j.

fft, with a single input argumen t x, computes the DFT of the input vector

or matrix. If rectangular array,

x is a vector, fft computes the DFT of the vector; if x is a

fft computes the DFT of each array column.

1-35

1 Filtering, Linear Systems and Transforms Overview

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]);

1-36

nd argument to

Aseco

esenting DFT length:

repr

y = fft(x,n);

is case,

In th

ncates the sequence if it is longer than

tru

he length of the input sequence. Execution time for

to t

gth,

len

cumentation for details about the algorithm.

fft pads the input sequence w ith zeros if it is shorter than n,or

n, of the DFT it performs; see the fft reference page in the MATLAB

fft specifies a number of points n for the transform,

n.Ifn isnotspecified,itdefaults

fft depends on the

Discrete Fourier Transform

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

x and 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, signal or image processing. The

fft2 and ifft2. These functions are useful for two-dimensional

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-37

1 Filtering, Linear Systems and Transforms Overview

1-38

Filter Design and Implementation

• “Filter Requirements and S pecification” 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

• “Selected Bibliography” on page 2-52

2 Filter 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 bee n 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 specif ied.

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, invokes an appropriate 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.

design method. To illustrate the object-oriented

fdesign,andthen

2-2

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

Filter Requirements and S pecification

[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, m ultiply 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

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

butter:

2-3

2 Filter Design and Implementation

IIR Filter Design

In this section...

“IIR vs. FIR Filters” on page 2-4

“Classical IIR F i lters ” 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 eliminates the nonlinear phase distortion of an IIR filter.

filtfilt function), which

2-4

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 class ical 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 s pe cification in terms of passband and stopband attenuation, and transition width(s).

IIR Filter Design

Other IIR Filter

The direct filte response approx multiband band

You can also us to design IIR f on page 7-15.

The generali section “Gen

r design function

ilters. These functions are discussed in “Parametric Modeling”

zed Butterworth design function

eralized Butterw orth Filter Design” on page 2-15.

IIR Filter M

The follow lists the f

Toolbox F

Filter Method Description Filter Functi ons

Analog Prototyping

ilters Methods and Available Functions

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.

ing table summarizes the various filter methods in the toolbox and

unctionsavailabletoimplementthesemethods.

yulewalk finds a filter with magnitude

imating a desired function. This is one way to create a

pass filter.

e the parametric modeling or syst em identification functions

maxflat is discussed in the

ethod Summary

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:

rect Des ign

Design digital filter directly in the discrete time-domain

lp2bp, lp2bs, lp2hp, lp2lp

Filter discretization functions:

bilinear, impinvar

yulewalk

2-5

2 Filter Design and Implementation

Toolbox Filters Methods and Available Functions (Continued)

Filter Method Description Filter Functi ons

by approximating a piecewise linear magnitude response.

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 ex ten s iv e 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.

2-6

Filter Design Functions

Filter Type Design Function

Bessel (analog only)

Butterworth

Chebyshev Type I

Chebyshev Type II

[b,a] = besself(n,Wn,options)

[z,p,k] = besself(n,Wn,options)

[A,B,C,D] = besself(n,Wn,options)

[b,a] = butter(n,Wn,options)

[z,p,k] = butter(n,Wn,options)

[A,B,C,D] = butter(n,Wn,options)

[b,a] = cheby1(n,Rp,Wn,options)

[z,p,k] = cheby1(n,Rp,Wn,options)

[A,B,C,D] = cheby1(n,Rp,Wn,options)

[b,a] = cheby2(n,Rs,Wn,options)

[z,p,k] = cheby2(n,Rs,Wn,options)

IIR Filter Design

[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 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 in normalized frequency (Nyquist

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

2-7

2 Filter Design and Implementation

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 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 SOS form with

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 d esign a wider passband and decimate the input signal.

zp2sos and then use the

dfilt to analyze or implement your filter.

2-8

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

Chebyshev Type I

Chebyshev Type II

Elliptic

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

[n,Wn] = buttord([1000 2000]/5000,[500 2500]/5000,1,60)

[n,Wn] = buttord(Wp,Ws,Rp,Rs)

[n,Wn] = cheb1ord(Wp, Ws, Rp, Rs)

[n,Wn] = cheb2ord(Wp, Ws, Rp, Rs)

[n,Wn] = ellipord(Wp, Ws, Rp, Rs)

butter function as follows.

IIR Filter Design

Wn =

0.1951 0.4080

[b,a] = butter(n,Wn);

An elliptic filter that meets the same requirements is given by

[n,Wn] = ellipord([1000 2000]/5000,[500 2500]/5000,1,60) n=

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; see Chapter 12, “Functions — Alphabetical List” for details.

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 any order N, the magnitude squared response has 2N-1 zero derivatives at these locations (maximally flat at

all, decreasing smoothly from

over

and ). Response is monotonic

to . at .

and ;for

2-9

2 Filter Design and Implementation

Chebyshev Type I Filter

The Cheby ideal and an equal r flat. The

Butterw

shev Type I filter minimizes th e absolute difference between the

actual frequency response over the entire passband by incorporating

ipple of Rp dB in the passband. Stopband response is maximally

transition from passband to stopband is more rapid than for the

orth filter.

at .

2-10

IIR Filter Design

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.

at .

2-11

2 Filter Design and Implementation

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 t ra nsition width.

at .

2-12

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 functio n .

Bessel filters generally require a higherfilterorderthanotherfiltersfor satisfactory stopband attenuation.

filter order n increases.

at and decreases as

IIR Filter Design

Note The lowpass filters shown above were created with the analog prototype functions find the zeros, poles, and gain of an order

besselap, buttap, cheb1ap, cheb2ap,andellipap. These functions

n analog 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 det ails.

To create similar plots, use n = 5 and, as needed, Rp = 0.5 and Rs = 20.For example, to create the elliptic filter plot:

k] = ellipap(5,0.5,20);

[z,p, w=lo h=fr semi

gspace(-1,1,1000); eqs(k*poly(z),poly(p),w);

logx(w,abs(h)), grid

2-13

2 Filter Design and Implementation

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 is intended specifically for filter design; “Parametric Modeling” on page 7-15 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. 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

yulewalk function, which

yulewalk’s name reflects its method for

2-14

[b,a] = yulewalk(n,f,m)

returns row vectors b and a containing the n+1 numerator and denominator coefficients of the order approximate those given in vectors ranging from 0 to 1, where 1 represents the Nyquist frequency. containing the desired magnitude response at the points in

n IIR filter whose frequency-magnitude characteristics

f and m. f is a vector of frequency points

m is a vector

f. f and 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 mag nitude response. See “FIR Filter Design” on page 2-17 for details.

Note that

yulewalk does not accept phase information, and no state m ents 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))

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 zero s. 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,

butter:

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

0.0317 0.0951 0.0951 0.0317

1.0000 -1.4590 0.9104 -0.1978

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

maxflat is just like the butter function, except

maxflat is the same as

2-15

2 Filter Design and Implementation

However, maxflat is more versatile because it allows you to design a filter with more zeros than poles:

[b,a] = maxflat(3,1,0.25) b=

The third input to maxflat is the half-power frequency, a frequency between 0 and 1 with a desired magnitude response of

You can also design linear phase filters that have the maximally flat property using the

0.0317 0.0951 0.0951 0.0317

1.0000 -1.4590 0.9104 -0.1978

0.0950 0.2849 0.2849 0.0950

1.0000 -0.2402

'sym' option:

2-16

maxflat(4,'sym',0.3) ans =

0.0331 0.2500 0.4337 0.2500 0.0331

For comple te details of the maxflat algorithm, see Selesnick and Burrus [2].

FIR Filter Design

In this section...

“FIR vs. IIR Filters” on pag e 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 Filter Design

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

2 Filter Design and Implementation

FIR Filter Summa

FIR Filters

Filter Design Method Description Filter Functions

Windowing

Multiband with Transi Bands

Constrain Least Squ

Arbitrary Response

Raised Cosine

Linea

pt for

Exce

ers only. The filter coefficients, or “taps,” of such filters obey either an even

filt

d symmetry relation. Depending on this symmetry, and on whether the

or od

r n of the filter is ev en or odd, a linear phase filter (stored in length n+1

orde

tor

vec

tion

ed ares

r Phase Filters

cfirpm, all of the FIR filter design functions design l inear phase

b) has certain inherent restrictions on its frequency response.

Apply window to truncated inverse Fourier transform of desired “brick wall” filter

Equiripple over sub-ba

Minimize squared integral error over entire frequency range subject to maximum error constraints

Arbitrary responses, including nonlinear phase and complex filters

Lowpass response with smooth, sinusoidal transition

or least squares approach

nds of the frequency range

fir1, fir2, kaiserord

firls, firpm, firpmord

fircls, f

cfirpm

firrcos

ircls1

2-18

FIR Filter Design

Linear Phase Filter Type

Filter Order S y mmetry of Coefficients

Type I Even

Type II

Odd

Type III Even

Type IV

Odd

Response

even:

Response H(f), f

= 0

No restriction

H(f), f (Nyquist)

No restriction

H(1)

= 1

= 0

restriction

odd:

H(0) = 0 H(1) = 0

H(0) = 0

No restriction

The phase delay and group delay of linear phase FIR filters are equal and constant over the frequency band. For an order n linear 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 design type I and II linear phase FIR filters by default. Both

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,

fir1, fir2, firls, firpm, fircls, fircls1,andfirrcos all

firls and

and nonlinear phase filters as well.

Note Because the frequency response of a type II filter is zero at the Nyquist frequency (“high” frequency), bandstop filters. For odd-valued

fir1 does not design type II highpass and

n in these cases, fir1 adds 1 to the order

and returns a type I filter.

2-19

2 Filter Design and Implementation

Windowing Metho

Consider the ide frequency of ω magnitude less between ω

This filter i noncausal. T applying a w this trunca filter wit

h a lowpass cutoff frequency ω

b = 0.4*sinc(0.4*(-25:25));

The windo theorem, filter, i the filt

fvtool(b,1)

this is the length 51 filter that bestapproximatestheideallowpass n the integrated least squares sense. The following command displays er’s frequency response in FVTool:

al, or “brick wall,” digital lowpass filter with a cutoff

ad/s. This filter has magnitude 1 at all frequencies with

than ω

. Its impulse response sequence h(n)is

and π

s not imp l em en table since its impulse response is infinite an d

o create a finite-duration impulse response, truncate it by

indow. By retaining the central section of impulse response in

tion, you obtain a linear phase FIR filter. For example, a length 51

w applied here is a simple rectangular window. By Parseval’s

, and magnitude 0 at frequencies with magnitude

of rad/s is

2-20

FIR Filter Design

Note that the y-axis shown in the f igure 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

2 Filter Design and Implementation

Note that the y-axis shown in the f igure below is in Magnitude Squared. You can set this by right-clicking on the axis label and selecting Magnitude

Squared from the menu.

2-22

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

fir1 and fir2 are based on this windowing process. Given a

FIR Filter Design

Standard Band F IR 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 b containing the coefficients of the order n Hamming-windowed filter. This is a lowpass, linear phase FIR filter with cutoff frequency the N yquist frequency, half the sam pling frequency. (Unlike other methods, here

Wn corresponds to the 6 dB point.) For a highpass filter, simply appen d

the string filter, specify frequencies; append the string

Wn. Wn is a number between 0 and 1, where 1 corresponds to

'high' to the function’s parameter l ist. For a bandpass or bandstop

Wn as a two-element vector containing the passband edge

'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,

Kaiser Window Order Estimation. The

fir1 applies a Hamming window.

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 w ell 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 respon se. This is in contrast to only designs filters in standard lowpas s, highpass, bandpass, and bandstop configurations.

The commands

fir1,which

2-23

2 Filter Design and Implementation

n=50; f = [0 .4 .5 1]; m = [1 1 0 0]; b = fir2(n,f,m);

return row vector b containing the n+1 coefficients of the order n FIR filter whose frequency-magnitude characteristics match those given by vectors and m. f is a vector of frequency points ranging from 0 to 1, where 1 represents the Nyquist frequency. response at the points specified in 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 design Hilbert transformers, diffe rentiators, 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 w hich the error is not minimized, and perform band dependent weighting of the minimization.

m is a v ector containing the desired magnitude

f. (The IIR counterpart of this function

fir1 and fir2 functions. These functions

2-24

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 appro ximation theory to design filters w ith 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.

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 am p litude 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 h av e an inherent trade-off similar to FIR design by windowing.

Tocompareleastsquarestoequiripplefilterdesign,use a similar filter. Type

bb = firls(n,f,a);

and compare their frequency respo n ses using FVTool:

fvtool(b,1,bb,1)

firls to create

2-25

2 Filter Design and Implementation

Note that the y-axis shown in the f igure below is in Magnitude Squared. You can set this by right-clicking on the axis label and selecting Magnitude

Squared from the menu.

2-26

lter designed with

The fi

rls

the but at the i erro poss

nk of frequency bands as lines over short frequency intervals.

Thi

fir

any and

f = [0 0.3 0.4 0.7 0.8 1]; % Band edges in pairs

filter has a better response over most of the passband and stopband,

the band edges (

deal than the

r over the passband and stopband is sm aller and, in fact, it is the smallest

ible for this band edge configuration and filter length.

use this scheme to represent any piecewise linear desired function with transition bands. bandstop filters; a bandpa ss example is

firpm exhibits equiripple behavior. Also note that

f = 0.4 and f = 0.5),theresponseisfurtherawayfrom

firpm filter. This shows that the firpm filter’s maximum

firpm and

firls and firpm design lowpass, highpass, bandpass,

FIR Filter Design

a = [0 0 1 1 0 0]; % Bandpass filter amplitude

Technically, these f and a vectors 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=[110 0 110 0 110 0 11];

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

2 Filter 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 f and a vectors; 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 H ilbert 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 x by 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 x is the complex signal that has x as its real part and the Hilbert transform of method (an alternative to the

hilbert function), you must delay x by half the

x as its imaginary part. For this FIR

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

2 Filter Design and Implementation

This method does not w ork directly for filters of odd order, which require a noninteger delay. In this case, the Transforms” on page 7-40, estimates the analytic signal. Alternatively, use the

resample function to dela y the signa l 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(ω) Approximate the ideal differentiator (with a delay) using a

'd' or 'differentiator' option:

b = firpm(21,[0 1],[0 pi],'d');

For a type III filter, t he 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');

hilbert function, described in “Specialized

= jω.

firpm or firls with

2-30

In the 'd' mode, firpm weights the error by 1/ω in nonzero amplitude bands to minimize the maximum relative error. nonzero amplitude bands in the

'd' mode.

firls weights the error by (1/ω)

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 s e veral 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

2 Filter 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 G ibbs’ phenomena.

There are two toolbox functions that implement this design technique.

Description Function

Constrained least square multiband FIR filter design

Constrained least square filter design for lowpass and highpass linear phase filters

fircls

fircls1

2-32

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 de sign 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.

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

2 Filter 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.4 9

• 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 respons e order 129 that meets these specifications:

2-34

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)

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 function enables you to specify the passband and stopband edges for the least squares weighting function, as well as a constant

k tha t 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:

fircls1

2-35

2 Filter 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

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)

fircls1,type

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

cfirpm a highly versatile and powerful

2-37

2 Filter 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 conve rgence to the optimal solution.

Multiband Filter Design

Consider a multiband filter with the following special frequency-domain characteristics.

Band Amplitude

[-1 -0.5] [5 1]

[-0.4 +0.3] [2 2]

[+0.4 +0.8] [2 1]

A linear-phase multiband filter may be designed using the predefined frequency-response function

Optimization Weighting

multiband,asfollows:

2-38

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

b = cfirpm(38,[-1 -0.5 -0.4 0.3 0.4 0.8], ...

[512221],[1105]);

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)

firpm:

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

2 Filter 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.

2-40

Filter Design with Reduced Delay

Consider the des ig n of a 62-t ap lowpass filter with a half-Nyquist cutoff. If we specify a negative offset value to 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)

lowpass filter design function, the

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

2 Filter 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.

2-42

If we c exact

29. U stop can a appl

ompare this nonlinear-phase filter to a linear-phase filter that has

ly 14.5 samples of group delay, the resulting filter is of order 2*14.5, or

sing

b = cfirpm(29,[0 0.5 0.55 1],'lowpass'), the passband and

band ripple is much g reate r for the order 29 filter. These comparisons

ssist you in deciding which filter is more appropriate for a specific

ication.

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 th e following steps.

1 Find an analog lowpass filter with cutoff frequency of 1 and translate this

prototype filter to the desired band configuration

2 Transform the filter to the digital domain.

Special Topics in IIR Filter Design

3 Discret

The too

ize the filter.

lbox provides functions for each of these steps.

Design Task Available functions

Analog lowpass

buttap, cheb1ap, besselap, ellipap, cheb2ap

prototype

ency

Frequ

sformation

tran

Discretization

lp2lp, lp2hp, lp2bp, lp2bs

near

bili

, impinvar

Alternatively, the butter, cheby1, cheb2ord, ellip,andbesself functions perform all steps of the filter design and the and

ellipord functions provide minimum order computation for IIR filters.

buttord, cheb1ord, cheb2ord,

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

2 Filter Design and Implementation

Analog Prototyp

This toolbox pro prototype filte approach to IIR

The table belo supported fil on page 2-4.

Filter Type Analog Prototype Function

Bessel

Butterworth

Chebyshev Type I

Chebyshev Type II

Elliptic

Freque

cond step in the analog prototyping design technique is the frequency

The se

formation of a lowpass prototype. The toolbox provides a set of functions

trans

nsform analog lowpass prototypes (with cutoff frequency of 1 rad/s)

to tra

andpass, highpass, bandstop, and lowpass filters of the desired cutoff

into b

uency.

freq

vides a number of functions to create lowpass analog

rs with cutoff frequency of 1, the firs t step in the classical

filter design.

w s ummarizes the analog prototype design functions for each

ter type; plots for each type are shown in “IIR Filter Design”

ncy Transformation

eDesign

[z,p,k] = besselap(n)

(n)

[z,p,k] = b

[z,p,k] = cheb1ap(n,Rp)

[z,p,k] = cheb2ap(n,Rs)

[z,p,k] = ellipap(n,Rp,Rs)

uttap

2-44

Frequency Transformation Transformation Function

Lowpass to lowpass

Lowpass to highpass

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

[At,Bt,Ct,Dt]

numt,dent]

[

At,Bt,Ct,Dt]

[

= lp2lp (A,B,C,D,Wo)

= lp2hp (num,den,Wo)

= lp2hp (A,B,C,D,Wo)

Special Topics in IIR Filter Design

Frequency Transformation Transformation Function

Lowpass to bandpass

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

= lp2bp (A,B,C,D,Wo,Bw)

= lp2bs( A,B,C,D,Wo,Bw)

Lowpass to ba

ndstop

[At,Bt,Ct,Dt]

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

[At,Bt,Ct,Dt]

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

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 substitution, so the output filter is twice the order of the input. For

lp2hp, the output filter is the same order as the input.

lp2bp and lp2bs, this is a second-order

lp2lp and

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(5,3);

Outputs z, p,andk contain the zeros, poles, and gain of a lowpass analog filter with cutoff frequen c y Ω

equal to 1 rad/s. Use the lp2bp function to

transform this lowpass prototype to a bandpass analog filter with band edges

2-45

2 Filter Design and Implementation

lp2bp function can accept it:

and . First, convert the filter to state-space form so the

[A,B,C,D] = zp2

ss(z,p,k); % Convert to state-space form.

Now, find the bandwidth and center frequency, and call lp2bp:

u1 = 0.1*2*pi; Bw = u2-u1; Wo = sqrt(u1* [At,Bt,Ct,D

u2 = 0.5*2*pi; % In radians per second

u2);

t] = lp2bp(A,B,C,D,Wo,Bw);

Finally, calculate the frequency response and plot its magnitude:

[b,a] = ss2t w = linspac h = freqs(b semilogy( xlabel('

f(At,Bt,Ct,Dt); % Convert to TF form e(0.01,1,500)*2*pi; % Generate frequency vector ,a,w); % Compute frequency response

w/2/pi,abs(h)), grid % Plot log magnitude vs. freq

Frequency (Hz)');

2-46

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

Mathworks SIGNAL PROCESSING TOOLBOX 6 user guide

Specifications and Main Features

Frequently Asked Questions

User Manual

Filtering, Linear Systems and Transforms Overview

Filter Implementation and Analysis

Filtering Overview

Convolution and Filtering

Filters and Transfer Functions

Filter Coefficients and Filter Names

Filtering with the filter Function

The filter Function

Other Functions for Filtering

Multirate Filter Bank Implementation

Anti-Causal, Zero-Phase Filter Implementation

Frequency Domain Filter Implementation

Impulse Response

Frequency Response

Digital Domain

Analog Domain

Delay

Zero-Pole Analysis

Linear System Models

Available Models

Discrete-Time System Models

Transfer Function

Zero-Pole-Gain

State-Space

Partial Fraction Expansion (Residue Form)

Second-Order Sections (SOS)

Lattice Structure

Convolution Matrix

Continuous-Time System Models

Linear System Transformations

Discrete Fourier Transform

Filter Design and Implementation

Filter Requirements and Specification

IIR Filter Design

IIRvs. FIRFilters

Classical IIR Filters

Classical IIR Filter Design Using Analog Prototyping

Complete Classical IIR Filter Design

Filter Design Functions

Designing IIR Filters to Frequency Domain Specifications

Comparison of Classical IIR Filter Types

Butterworth Filter

Chebyshev Type I Filter

Chebyshev Type II Filter

Elliptic Filter

Bessel Filter

Direct IIR Filter Design

Generalized Butterworth Filter Design

FIR Filter Design

FIRvs. IIRFilters

FIR Filters

Standard Band F IR Filter Design: fir1

Multiband FIR Filter Design: fir2

Multiband FIR Filter Design with Transition Bands

Basic Configurations

The Weight Vector

Anti-Symmetric Filters / Hilbert Transformers

Differentiators

Constrained Least Squares FIR Filter Design

Basic Lowpass and Highpass CLS Filter Design

Multiband CLS Filter Design

Weighted CLS Filter Design

Arbitrary-Response Filter Design

Multiband Filter Design

Filter Design with Reduced Delay

Special Topics in IIR Filter Design

Classic IIR Filter Design

Frequency Transformation Transformation Function

Filter Discretization