《DSP using MATLAB》Problem 7.6

代码：

子函数ampl_res

function [Hr,w,P,L] = ampl_res(h);
%
% function [Hr,w,P,L] = Ampl_res(h)
% Computes Amplitude response Hr(w) and its polynomial P of order L,
% given a linear-phase FIR filter impulse response h.
% The type of filter is determined automatically by the suroutine.
%
% Hr = Amplitude Response
% w = Frequence between [0 pi] over which Hr is computed
% P = Polynomial coefficients
% L = Order of P
% h = Linear Phase filter impulse response
%
M = length(h);
if rem(M,2)==1
   if all(h(1:(M-1)/2)==h(M:-1:(M+3)/2)) [Hr,w,P,L] = Hr_Type1(h);
   elseif all(h(1:(M-1)/2)==-h(M:-1:(M+3)/2)) & h((M+1)/2)==0, [Hr,w,P,L] = Hr_Type3(h); 
       else disp('not a linear-phase filter, check h'), return, 
   end
elseif all(h(1:M/2)==h(M:-1:M/2+1)); [Hr,w,P,L] = Hr_Type2(h);
   elseif all(h(1:M/2)==-h(M:-1:M/2+1)); [Hr,w,P,L] = Hr_Type4(h);
       else disp('not a linear-phase filter, check h'), return, 
   end
end

　　第1小题代码

%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%%            Output Info about this m-file
fprintf('
***********************************************************
');
fprintf('        <DSP using MATLAB> Problem 7.6 

');

banner();
%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

n1 = 0; n2 = 11;
n = [n1:n2];
h1n = 0.9 .^ abs(n-5) .* cos(pi*(n-5)/12) .* ( stepseq(0,n1,n2) - stepseq(11,n1,n2) );

h1n = h1n(1:11);

figure('NumberTitle', 'off', 'Name', 'Problem 7.6')
set(gcf,'Color','white'); 
stem([0:10], h1n); grid on;
xlabel('n'); ylabel('h1(n)'); title('Impulse Response');

[db, mag, pha, grd, w] = freqz_m(h1n, 1);
[Hr, ww, a, P] = ampl_res(h1n);
%[Hr, ww, a, P] = Hr_Type1(h1n);


% Plot 

figure('NumberTitle', 'off', 'Name', 'Problem 7.6')
set(gcf,'Color','white'); 

subplot(2,2,1); plot(w/pi, db); grid on; axis([0 2 -100 10]); 
xlabel('frequency in pi units'); ylabel('Decibels'); title('Magnitude Response in dB');
set(gca,'YTickMode','manual','YTick',[-80,-40,0]);
set(gca,'YTickLabelMode','manual','YTickLabels',['80';'40';' 0']);

subplot(2,2,3); plot(w/pi, mag); grid on; %axis([0 1 -100 10]); 
xlabel('frequency in pi units'); ylabel('Absolute'); title('Magnitude Response in absolute');
subplot(2,2,2); plot(w/pi, pha); grid on; %axis([0 1 -100 10]); 
xlabel('frequency in pi units'); ylabel('Rad'); title('Phase Response in Radians');
subplot(2,2,4); plot(w/pi, grd*pi/180);  grid on; %axis([0 1 -100 10]); 
xlabel('frequency in pi units'); ylabel('Rad'); title('Group Delay');


figure('NumberTitle', 'off', 'Name', 'Problem 7.6')
set(gcf,'Color','white'); 

subplot(2,1,1); plot(ww/pi, Hr); grid on; axis([0, 1, 0, 8]); 
xlabel('frequency in pi nuits'); ylabel('Hr'); title('Amplitude Response');

subplot(2,1,2); stem([n1:n2-1], h1n); grid on;%axis([-1, M, -1.1, 1.1]);
xlabel('n'); ylabel('h1(n)'); title('Impulse Response');

　　运行结果：