《DSP using MATLAB》Problem 7.37

代码：

%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

%%            Output Info about this m-file

fprintf('\n***********************************************************\n');

fprintf('        <DSP using MATLAB> 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

%weights

fprintf('\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 order

M = N+1

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;

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方法要比窗函数法、频率采样法设计的滤波器长度小，节约硬件资源，这里没有证明，以后有时间再算吧。