代码：
%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%% Output Info about this m-filefprintf('\n***********************************************************\n');fprintf(' Problem 7.37 \n\n');banner();%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++w1 = 0; w2 = 0.30*pi; delta1 = 0.05; gain1 = 0.95;w3 = 0.35*pi; w4 = 0.50*pi; delta2 = 0.10; gain2 = 2.00;w5 = 0.55*pi; w6 = 0.75*pi; delta3 = 0.05; gain3 = 3.05;w7 = 0.80*pi; w8 = pi; delta4 = 0.05; gain4 = 4.05;f = [w2, w3, w4, w5, w6, w7]/pi; m = [1, 1, 1, 1]; delta = [delta1, delta2, delta3, delta4];[N, f0, m0, weights] = firpmord(f, m, delta); fprintf('\n-------------------------------------------\n');fprintf('\n Results by firpmord function:\n');N%f0%m0%weightsfprintf('\n-------------------------------------------\n');weights = [delta2/delta1 1 delta2/delta3 delta2/delta4]deltaH = max([delta1,delta2,delta3,delta4]); deltaL = min([delta1,delta2,delta3,delta4]);%Dw1 = min((w3-w2), (w5-w4));%Dw2 = min((w5-w4), (w7-w6));%Dw = min(Dw1, Dw2);%M = ceil((-20*log10((delta1*delta2*delta3*delta4)^(1/4)) - 13) / (2.285*Dw) + 1)f = [ 0 w2 w3 w4 w5 w6 w7 pi]/pi;m = [ 0.95 0.95 2 2 3.05 3.05 4.05 4.05];h = firpm(N, f, m, weights); % even-number orderM = N+1delta_w = 2*pi/1000;[db, mag, pha, grd, w] = freqz_m(h, [1]);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;w7i = floor(w7/delta_w)+1; w8i = floor(w8/delta_w)+1;Asd1 = -max(db(w1i:w2i))%Asd2 = -max(db(w3i:w4i))%Asd3 = -max(db(w5i:w6i))%Asd4 = -max(db(w7i:w8i))%[Hr, ww, a, L] = Hr_Type1(h);[Hr,omega,P,L] = ampl_res(h);l = 0:M-1;%% -------------------------------------------------%% Plot%% ------------------------------------------------- figure('NumberTitle', 'off', 'Name', 'Problem 7.37')set(gcf,'Color','white'); subplot(2,2,1); stem(l, h); axis([-1, M, -1.0, 2.4]); grid on;xlabel('n'); ylabel('h(n)'); title('Actual Impulse Response, M=45');set(gca,'XTickMode','manual','XTick',[0,(M-1)/2,M-1])set(gca,'YTickMode','manual','YTick',[-1.0:0.4:2.4])subplot(2,2,2); plot(w/pi, db); axis([0, 1, -20, 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',[-20,-12,-5,-2,0]);set(gca,'YTickLabelMode','manual','YTickLabel',['20';'12';' 5';' 2';' 0']);subplot(2,2,3); plot(omega/pi, Hr); axis([0, 1, -0.2, 4.5]); grid on;xlabel('frequency in \pi nuits'); ylabel('Hr(\omega)'); title('Amplitude Response');set(gca,'XTickMode','manual','XTick',f)set(gca,'YTickMode','manual','YTick',[0.9,1.0,1.9,2.1,3.0,3.1,4.0,4.1]);delta_w = 2*pi/1000;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);b4w = omega(w7i:w8i)/pi; b4e = (Hr(w7i:w8i)-m(7)); %b4e = (Hr(w7i:w8i)-m(7))*weights(4);plot(b1w,b1e,b2w,b2e,b3w,b3e,b4w,b4e); 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.37 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(\omega)'); title('Amplitude Response');set(gca,'YTickMode','manual','YTick',[0.9,1.0,1.9,2.1,3.0,3.1,4.0,4.1]);set(gca,'XTickMode','manual','XTick',f); 运行结果：
先用firpmord求出N、f、m、weights，这些参数作为firpm函数参数，再求出脉冲响应序列
求出脉冲响应序列的幅度谱、误差响应函数
振幅谱，按照题目要求画了网格线
依前面经验，Parks-McClellan方法要比窗函数法、频率采样法设计的滤波器长度小，节约硬件资源，这里没有证明，以后有时间再算吧。