Sign in

username:

password:



Not a member?

Search Online Books



Search tips

Free Online Books

Sponsor

Industry's highest performing at the lowest power DSPs now as low as $5.00*
Start development today!
*volume pricing for 10ku

Chapters

See Also

Embedded SystemsFPGAElectronics

Introduction to Digital Filters
    Frequency Response Analysis
       Frequency Response as a Ratio of DTFTs
          Frequency Response in Matlab

Chapter Contents:

Search Introduction to Digital Filters

  

Book Index | Global Index

Would you like to be notified by email when Julius Orion Smith III publishes a new entry into his blog?

  

Frequency Response in Matlab

In practice, we usually work with a sampled frequency axis. That is, instead of evaluating the transfer function $ H(z)=B(z)/A(z)$ at $ z=e^{j\omega T}$ to obtain the frequency response $ H(e^{j\omega T})$, where $ \omega$ is continuous radian frequency, we compute instead

$\displaystyle H(e^{j\omega_k T}) \eqsp \frac{B(e^{j\omega_k T})}{A(e^{j\omega_k T})},\quad e^{j\omega_k T}\isdefs e^{j2\pi k / N_s},\quad k=0,1,2,\ldots,N_s-1, $

where $ N_s$ is the desired number of spectral samples around the unit circle in the $ z$ plane. From Eq.$ \,$(7.5) we have that this is the same thing as the sampled DTFT of $ {\underline{b}}$ divided by the sampled DTFT of $ {\underline{a}}$:

$\displaystyle H(e^{j\omega_k T}) \eqsp \frac{\mbox{{\sc DTFT}}_{\omega_k T}({\underline{b}})}{\mbox{{\sc DTFT}}_{\omega_k T}({\underline{a}})} $

The uniformly sampled DTFT has its own name: the discrete Fourier transform (DFT) [84]. Thus, we can write

$\displaystyle H(e^{j\omega_k T}) \eqsp \frac{\mbox{{\sc DFT}}_{\omega_k T}({\un... ...}{\mbox{{\sc DFT}}_{\omega_k T}({\underline{a}})}, \quad k=0,1,2,\ldots,N_s-1 $

where

   DFT$\displaystyle _{\omega_k T}(x) \isdefs \sum_{n=0}^{N_s-1} x(n) e^{-j\omega_k nT} $

and $ \omega_k \isdef 2\pi f_sk/N_s$, where $ f_s=1/T$ denotes the sampling rate in Hz.

To avoid undersampling $ B(e^{j\omega T})$, we must have $ N_s\ge M$, and to avoid undersampling $ A(e^{j\omega T})$, we must have $ N_s\ge N$. In general, $ H(e^{j\omega T})$ will be undersampled (when $ N>0$), because it is the quotient of $ B(e^{j\omega T})$ over $ A(e^{j\omega T})$. This means, for example, that computing the impulse response $ h(n)$ from the sampled frequency response $ H(e^{j\omega_k T})$ will be time aliased in general. I.e.,

$\displaystyle h(n) \eqsp$   IDFT$\displaystyle _n(H) \isdefs \frac{1}{N_s}\sum_{k=0}^{N_s-1} H(e^{j\omega_k T})e^{j\omega_k nT} $

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 $ N_s$ 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 $ N_s>7/(1-R_{\mbox{max}})$, where $ R_{\mbox{max}}$ denotes the maximum pole magnitude. This choice provides more than 60 dB of decay in the impulse response over a duration of $ N_s$ samples, which is the time-aliasing block size. (The time to decay 60 dB, or ``$ t_{60}$'', is a little less than $ 7$ time constants [84], and the time-constant of decay for a single pole at radius $ R$ can be approximated by $ 1/(1-R)$ samples, when $ R$ is close to 1, as derived in §8.6.)

As is well known, when the DFT length $ N_s$ is a power of 2, e.g., $ N_s=2^{10}=1024$, 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 $ N_s$ is a power of 2).

Figure 7.1: Matlab function for computing and optionally plotting the frequency response of an IIR digital filter.
  
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

Previous: Frequency Response as a Ratio of DTFTs
Next: Example LPF Frequency Response Using freqz

Order a Hardcopy of Introduction to Digital Filters

About the Author: Julius Orion Smith III
Julius Smith's background is in electrical engineering (BS Rice 1975, PhD Stanford 1983). He is presently Professor of Music and Associate Professor (by courtesy) of Electrical Engineering at Stanford's Center for Computer Research in Music and Acoustics (CCRMA), teaching courses and pursuing research related to signal processing applied to music and audio systems. See http://ccrma.stanford.edu/~jos/ for details.


Comments

No comments yet for this page


Add a Comment
You need to login before you can post a comment (best way to prevent spam). ( Not a member? )