Practical Frequency-Response Analysis

The preceding examples were constructed to be tutorial on the level of this book. However, to provide a more practical perspective, a more typical frequency response analysis is presented in Fig.2.11.

A comparison of computed and theoretical frequency response curves is shown in Fig.2.12. There is no visible difference, and the only source of error is computational round-off error. Not only is this method as accurate as any other, it is by far the fastest, because it uses the Fast Fourier Transform (FFT).

This FFT method for computing the frequency response is based on the fact that the frequency response equals the filter transfer function $H(z)=B(z)/A(z)$ evaluated on the unit circle $z=e^{j\omega T}$ in the complex $z$ plane. We will get to these concepts later. For now, just note the ease with which we can compute the frequency response numerically in matlab. In fact, the length $N$ frequency response of the simplest lowpass filter $y(n) = x(n) + x(n - 1)$ can be computed using a single line of matlab code:

        H = fft([1,1],N);

When $N$ is a power of 2, the radix-2 FFT algorithm can be used for high-speed execution.^3.6 When $N>2$, the FFT input is automatically ``zero padded'' in the time domain, resulting in interpolation of the frequency response [83].^3.7 In other words, the above line of code is equivalent to

        H = fft([1,1,zeros(1,N-2)]);

when $N\ge 2$. The code in Fig.2.11 carries out explicit zero-padding for clarity.

In both Matlab and Octave, there is a built-in function freqz which uses this FFT method for calculating the frequency response for almost any digital filter $B(z)/A(z)$ (any causal, finite-order, linear, and time-invariant digital filter, as we will later see).

Figure: Main program in matlab for finding the frequency response of the simplest low-pass filter $y(n) = x(n) + x(n - 1)$ using the FFT. See §H.13 for the final plotting utility.

 

% simplpnfa.m - matlab program for frequency analysis 
%               of the simplest lowpass filter:
%
%                    y(n) = x(n)+x(n-1)}
%
%               the way people do it in practice.

B = [1,1]; % filter feedforward coefficients 
A = 1;     % filter feedback coefficients

N=128;     % FFT size = number of COMPLEX sinusoids
fs = 1;    % sampling rate in Hz (arbitrary)

Bzp = [B, zeros(1,N-length(B))]; % zero-pad for the FFT

H = fft(Bzp);   % length(Bzp) should be a power of 2

if length(A)>1  % we're not using this here,
  Azp = [A,zeros(1,N-length(A))]; % but show it anyway.
  % [Should guard against fft(Azp)==0 for some index]
  H = H ./ fft(A,N); % denominator from feedback coeffs
end

% Discard the frequency-response samples at 
% negative frequencies, but keep the samples at
% dc and fs/2:

nnfi = (1:N/2+1);     % nonnegative-frequency indices
Hnnf = H(nnfi);       % lose negative-frequency samples
nnfb = nnfi-1;        % corresponding bin numbers
f = nnfb*fs/N;        % frequency axis in Hz
gains = abs(Hnnf);    % amplitude response
phases = angle(Hnnf); % phase response

plotfile = '../eps/simplpnfa.eps';
swanalmainplot;    % final plots and comparison to theory

Figure 2.12: Overlay of theoretical and numerically determined frequency response using the FFT method.