《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('\n***********************************************************\n');

fprintf('        <DSP using MATLAB> Problem 7.6 \n\n');

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');