代码：
%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%% Output Info about this m-file
fprintf('
***********************************************************
');
fprintf(' <DSP using MATLAB> Problem 7.35
');
banner();
%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
w1 = 0; w2 = 0.20*pi; delta1 = 0.20; gain1 = 2.0;
w3 = 0.25*pi; w4 = 0.45*pi; delta2 = 0.05; gain2 = 0.0;
w5 = 0.55*pi; w6 = 0.70*pi; delta3 = 0.30; gain3 = 3.0;
w7 = 0.75*pi; w8 = pi; delta4 = 0.10; gain4 = 1.0;
f = [w2, w3, w4, w5, w6, w7]/pi; m = [1, 0, 1, 1]; delta = [delta1, delta2, delta3, delta4];
[N, f0, m0, weights] = firpmord(f, m, delta);
fprintf('
-------------------------------------------
');
fprintf('
Results will be used by firpm function:
');
N
%f0
%m0
%weights
fprintf('
-------------------------------------------
');
weights = [delta3/delta1 delta3/delta2 1 delta3/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 = [ 2 2 0 0 3 3 1 1];
h = firpm(N+4, f, m, weights); % even-number order
M = N+5
delta_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;
Asd = -max(db(w3i:w4i))
%[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.35')
set(gcf,'Color','white');
subplot(2,2,1); stem(l, h); axis([-1, M, -0.45, 1.45]); grid on;
xlabel('n'); ylabel('h(n)'); title('Actual Impulse Response, M=35');
set(gca,'XTickMode','manual','XTick',[0,M-1])
set(gca,'YTickMode','manual','YTick',[-0.4:0.2:1.4])
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,-34,-30,0]);
set(gca,'YTickLabelMode','manual','YTickLabel',['60';'34';'30';' 0']);
subplot(2,2,3); plot(omega/pi, Hr); axis([0, 1, -0.2, 4.0]); 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.05,0.05,0.9,1.1,1.8,2.2,2.7,3.3]);
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.35 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(Not to Scale)'); title('Amplitude Response');
set(gca,'YTickMode','manual','YTick',[-0.05,0.05,0.9,1.1,1.8,2.2,2.7,3.3]);
%set(gca,'YTickLabelMode','manual','YTickLabel',['90';'40';' 0']);
set(gca,'XTickMode','manual','XTick',[0,0.2,0.25,0.45,0.55,0.7,0.75,1]);
运行结果：
30阶时，稍有欠缺
增大滤波器长度，M=35时，
幅度谱、误差函数如下
振幅谱如下
