DSP using MATLAB示例Example3.16

代码：

b = [0.0181,  0.0543, 0.0543,  0.0181];         % filter coefficient array b
a = [1.0000, -1.7600, 1.1829, -0.2781];         % filter coefficient array a
m = 0:length(b)-1; l = 0:length(a)-1;           % index arrays m and l
K = 500; k = 0:1:K;                             % index array k for frequencies
w = k*pi/500;                                   % freqency between 0 and +pi, [0,pi] axis divided into 501 points.

num = b * exp(-j*m'*w);                         % Numerator calculations
den = a * exp(-j*l'*w);                         % Denominator claculations
H = num ./ den;                                 % Frequency response

magH  = abs(H);  angH  = angle(H);  realH  = real(H);  imagH  = imag(H);

%% --------------------------------------------------------------------
%%              START H's  mag ang real imag
%% --------------------------------------------------------------------
figure('NumberTitle', 'off', 'Name', 'Example3.16');
set(gcf,'Color','white'); 
subplot(2,2,1); plot(w/pi,magH); grid on;  axis([0,1,0,1.5]); 
title('Magnitude Response');
xlabel('frequency in pi units'); ylabel('Magnitude  |H|'); 
subplot(2,2,3); plot(w/pi, angH/pi); grid on; % axis([-1,1,-1,1]);
title('Phase Response');
xlabel('frequency in pi units'); ylabel('Radians/pi');

subplot('2,2,2'); plot(w/pi, realH); grid on;
title('Real Part');
xlabel('frequency in pi units'); ylabel('Real');
subplot('2,2,4'); plot(w/pi, imagH); grid on;
title('Imaginary Part');
xlabel('frequency in pi units'); ylabel('Imaginary');
%% -------------------------------------------------------------------
%%             END X's  mag ang real imag
%% -------------------------------------------------------------------

运行结果：

　　

    上图看出，是个低通滤波器。