《DSP using MATLAB》Problem 7.32

代码：

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

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

w1 = 0;      w2 = 0.4*pi; delta1 = 0.05;  gain1 = 0.4;
w3 = 0.5*pi; w4 = 0.7*pi; delta2 = 0.05;  gain2 = 0.95;
w5 = 0.8*pi; w6 = pi;     delta3 = 0.025; gain3 = 0.0;

f = [w2, w3, w4, w5]/pi; m = [1, 1, 0]; delta = [delta1, delta2, delta3];

[N, f0, m0, weights] = firpmord(f, m, delta); 
fprintf('\n-------------------------------------------\n');
fprintf('\n Results by firpmord function:\n');
N
%f0
%m0
%weights
fprintf('\n-------------------------------------------\n');


weights = [1  delta1/delta2 delta1/delta3]
deltaH = max([delta1,delta2,delta3]); deltaL = min([delta1,delta2,delta3]);

%Dw = min((w3-w2), (w5-w3));
%M = ceil((-20*log10((delta1*delta2*delta3)^(1/3)) - 13) / (2.285*Dw) + 1)

f = [ 0   w2    w3    w4    w5     pi]/pi;
m = [ 0.4 0.4   0.95  0.95  0.025  0.025];

h = firpm(N+1, f, m, weights);               % even-number order
[db, mag, pha, grd, w] = freqz_m(h, [1]);
delta_w = 2*pi/1000; 
w1i = floor(w1/delta_w)+1; w2i = floor(w2/delta_w)+1;
w3i = floor(w3/delta_w)+1; w4i = floor(w4/delta_w)+1;
w5i = floor(w5/delta_w)+1; w6i = floor(w6/delta_w)+1;

Asd = -max(db(w5i:w6i))

%[Hr, ww, a, L] = Hr_Type1(h);
[Hr,omega,P,L] = ampl_res(h);

M = N+2
l = 0:M-1;

%% -------------------------------------------------
%%                    Plot
%% -------------------------------------------------  

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

subplot(2,2,1); stem(l, h); axis([-1, M, -0.3, 0.6]); grid on;
xlabel('n'); ylabel('h(n)'); title('Actual Impulse Response, M=27');
set(gca,'XTickMode','manual','XTick',[0,M-1])
set(gca,'YTickMode','manual','YTick',[-0.3:0.1:0.6])

subplot(2,2,2); plot(w/pi, db); axis([0, 1, -80, 10]); grid on;
xlabel('frequency in \pi units'); ylabel('Decibels'); title('Magnitude Response in dB ');
set(gca,'XTickMode','manual','XTick',f)
set(gca,'YTickMode','manual','YTick',[-60,-26,0]);
set(gca,'YTickLabelMode','manual','YTickLabel',['60';'26';' 0']);

subplot(2,2,3); plot(omega/pi, Hr); axis([0, 1, -0.2, 1.2]); grid on;
xlabel('frequency in \pi nuits'); ylabel('Hr(w)'); title('Amplitude Response');
set(gca,'XTickMode','manual','XTick',f)
set(gca,'YTickMode','manual','YTick',[0,0.05,0.35,0.45,0.9,1]);

subplot(2,2,4); axis([0, 1, -deltaH, deltaH]); 
b1w = omega(w1i:w2i)/pi; b1e = (Hr(w1i:w2i)-m(1)); %b1e = (Hr(w1i:w2i)-m(1))*weights(1);
b2w = omega(w3i:w4i)/pi; b2e = (Hr(w3i:w4i)-m(3)); %b2e = (Hr(w3i:w4i)-m(3))*weights(2);
b3w = omega(w5i:w6i)/pi; b3e = (Hr(w5i:w6i)-m(5)); %b3e = (Hr(w5i:w6i)-m(5))*weights(3);
plot(b1w,b1e,b2w,b2e,b3w,b3e); grid on;
xlabel('frequency in \pi units'); ylabel('Hr(w)'); title('Error Response'); %title('Weighted Error');
set(gca,'XTickMode','manual','XTick',f);
set(gca,'YTickMode','manual','YTick',[-deltaH,0,deltaH]);


figure('NumberTitle', 'off', 'Name', 'Problem 7.32 AmpRes of h(n), Parks-McClellan Method')
set(gcf,'Color','white'); 

plot(omega/pi, Hr); grid on; %axis([0 1 -100 10]); 
xlabel('frequency in \pi units'); ylabel('Hr'); title('Amplitude Response');
set(gca,'YTickMode','manual','YTick',[0,delta2, 0.35, 0.45, 0.95-delta1, 0.95, 0.95+delta1]);
%set(gca,'YTickLabelMode','manual','YTickLabel',['90';'40';' 0']);
set(gca,'XTickMode','manual','XTick',[0,0.4,0.5,0.7,0.8,1]);

　　运行结果：