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 evaluated on the unit circle in the complex 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 frequency response of the simplest lowpass filter can be computed using a single line of matlab code:
H = fft([1,1],N);When is a power of 2, the radix-2 FFT algorithm can be used for high-speed execution.3.6 When , the FFT input is automatically ``zero padded'' in the time domain, resulting in interpolation of the frequency response .3.7 In other words, the above line of code is equivalent to
H = fft([1,1,zeros(1,N-2)]);when . 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 (any causal, finite-order, linear, and time-invariant digital filter, as we will later see).
% 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