## Frequency Response as a Ratio of DTFTs

From Eq.(6.5), we have , so that the frequency response isThis expression provides a convenient basis for the computation of an IIR frequency response in software, as we pursue further in the next section.

### Frequency Response in Matlab

In practice, we usually work with a*sampled*frequency axis. That is, instead of evaluating the transfer function at to obtain the frequency response , where is

*continuous*radian frequency, we compute instead

*sampled DTFT*of divided by the

*sampled DTFT*of :

*discrete Fourier transform (DFT)*[84]. Thus, we can write

DFT

and
, where denotes the
sampling rate in Hz.
To avoid undersampling
, we must have , and to
avoid undersampling
, we must have . In general,
*will be undersampled*(when ), because it is the quotient of over . This means, for example, that computing the impulse response from the sampled frequency response will be

*time aliased*in general.

*I.e.*,

IDFT

will be time-aliased in the IIR case. In other words, an infinitely
long impulse response cannot be Fourier transformed using a
finite-length DFT, and this corresponds to not being able to sample
the frequency response of an IIR filter without some loss of
information. In practice, we simply choose sufficiently large
so that the sampled frequency response is accurate enough for our
needs. A conservative practical rule of thumb when analyzing stable
digital filters is to choose
, where
denotes the maximum pole magnitude. This choice
provides more than 60 dB of decay in the impulse response over a
duration of samples, which is the time-aliasing block size.
(The time to decay 60 dB, or ``'', is a little less than
time constants [84], and the time-constant of decay for a
single pole at radius can be approximated by samples,
when is close to 1, as derived in §8.6.)
As is well known, when the DFT length is a power of 2, *e.g.*, , the DFT can be computed extremely efficiently using the

*Fast Fourier Transform (FFT)*. Figure 7.1 gives an example matlab script for computing the frequency response of an IIR digital filter using two FFTs. The Matlab function

`freqz`also uses this method when possible (

*e.g.*, when is a power of 2).

function [H,w] = myfreqz(B,A,N,whole,fs) %MYFREQZ Frequency response of IIR filter B(z)/A(z). % N = number of uniform frequency-samples desired % H = returned frequency-response samples (length N) % w = frequency axis for H (length N) in radians/sec % Compatible with simple usages of FREQZ in Matlab. % FREQZ(B,A,N,whole) uses N points around the whole % unit circle, where 'whole' is any nonzero value. % If whole=0, points go from theta=0 to pi*(N-1)/N. % FREQZ(B,A,N,whole,fs) sets the assumed sampling % rate to fs Hz instead of the default value of 1. % If there are no output arguments, the amplitude and % phase responses are displayed. Poles cannot be % on the unit circle. A = A(:).'; na = length(A); % normalize to row vectors B = B(:).'; nb = length(B); if nargin < 3, N = 1024; end if nargin < 4, if isreal(b) & isreal(a), whole=0; else whole=1; end; end if nargin < 5, fs = 1; end Nf = 2*N; if whole, Nf = N; end w = (2*pi*fs*(0:Nf-1)/Nf)'; H = fft([B zeros(1,Nf-nb)]) ./ fft([A zeros(1,Nf-na)]); if whole==0, w = w(1:N); H = H(1:N); end if nargout==0 % Display frequency response if fs==1, flab = 'Frequency (cyles/sample)'; else, flab = 'Frequency (Hz)'; end subplot(2,1,1); % In octave, labels go before plot: plot([0:N-1]*fs/N,20*log10(abs(H)),'-k'); grid('on'); xlabel(flab'); ylabel('Magnitude (dB)'); subplot(2,1,2); plot([0:N-1]*fs/N,angle(H),'-k'); grid('on'); xlabel(flab); ylabel('Phase'); end |

###
Example LPF Frequency Response Using `freqz`

Figure 7.2 lists a short matlab program illustrating usage of
`freqz`in Octave (as found in the

`octave-forge`package). The same code should also run in Matlab, provided the Signal Processing Toolbox is available. The lines of code not pertaining to plots are the following:

[B,A] = ellip(4,1,20,0.5); % Design lowpass filter B(z)/A(z) [H,w] = freqz(B,A); % Compute frequency response H(w)The filter example is a recursive fourth-order elliptic function lowpass filter cutting off at half the Nyquist limit (``'' in the fourth argument to

`ellip`). The maximum passband ripple

^{8.2}is 1 dB (2nd argument), and the maximum stopband ripple is 20 dB (3rd arg). The sampled frequency response is returned in the

`H`array, and the specific radian frequency samples corresponding to

`H`are returned in the

`w`(``omega'') array. An immediate plot can be obtained in either Matlab or Octave by simply typing

plot(w,abs(H)); plot(w,angle(H));However, the example of Fig.7.2 uses more detailed ``compatibility'' functions listed in Appendix J. In particular, the

`freqplot`utility is a simple compatibility wrapper for

`plot`with label and title support (see §J.2 for Octave and Matlab version listings), and

`saveplot`is a trivial compatibility wrapper for the

`freqplot`plots are shown in Fig.7.3(a) and Fig.7.3(b).

^{8.3}

[B,A] = ellip(4,1,20,0.5); % Design the lowpass filter [H,w] = freqz(B,A); % Compute its frequency response % Plot the frequency response H(w): % figure(1); freqplot(w,abs(H),'-k','Amplitude Response',... 'Frequency (rad/sample)', 'Gain'); saveplot('../eps/freqzdemo1.eps'); figure(2); freqplot(w,angle(H),'-k','Phase Response',... 'Frequency (rad/sample)', 'Phase (rad)'); saveplot('../eps/freqzdemo2.eps'); % Plot frequency response in a "multiplot" like Matlab uses: % figure(3); plotfr(H,w/(2*pi)); if exist('OCTAVE_VERSION') disp('Cannot save multiplots to disk in Octave') else saveplot('../eps/freqzdemo3.eps'); end |

Amplitude
Response
Phase Response |

**Next Section:**

Phase and Group Delay

**Previous Section:**

Polar Form of the Frequency Response