《DSP using MATLAB》Problem 7.29

代码：

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

');

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

% bandpass
ws1 = 0.25*pi; wp1 = 0.4*pi; wp2 = 0.6*pi; ws2 = 0.75*pi;
As = 50; Rp = 1;

[delta1, delta2] = db2delta(Rp, As);
deltaH = max(delta1,delta2); deltaL = min(delta1,delta2);

f = [ws1, wp1, wp2, ws2]/pi; m = [0, 1, 0]; delta = [delta2, delta1, delta2];

[N, f, m, weights] = firpmord(f, m, delta); 
N

h = firpm(N, f, m, weights);
[db, mag, pha, grd, w] = freqz_m(h, [1]);
delta_w = 2*pi/1000; 
ws1i = floor(ws1/delta_w)+1; wp1i = floor(wp1/delta_w)+1;
ws2i = floor(ws2/delta_w)+1; wp2i = floor(wp2/delta_w)+1;

Asd = -max(db(1:1:ws1i))

N = N + 2            % Length M=25
h = firpm(N, f, m, weights);

[db, mag, pha, grd, w] = freqz_m(h, [1]);
Asd = -max(db(1:1:ws1i))

M = N + 1
l = 0:M-1;
%% --------------------------------------------------
%%                   Type-1 BPF
%% --------------------------------------------------  
[Hr, ww, a, L] = Hr_Type1(h);

Rp = -(min(db(floor(wp1/delta_w)+1 :1: floor(wp2/delta_w))));     % Actual Passband Ripple
fprintf('
Actual Passband Ripple is %.4f dB.
', Rp);

As = -round(max(db(floor(ws2/delta_w)+1 : 1 : 501)));        % Min Stopband attenuation
fprintf('
Min Stopband attenuation is %.4f dB.
', As);

[delta1_db, delta2_db] = db2delta(Rp, As)


% Plot
figure('NumberTitle', 'off', 'Name', 'Problem 7.29 h(n), Parks-McClellan Method')
set(gcf,'Color','white'); 
subplot(2,2,1); stem([0:M-1], h); axis([0 M-1 -0.5 0.5]); grid on;
xlabel('n'); ylabel('h(n)'); title('Actual Impulse Response, M=25');

subplot(2,2,2); plot(w/pi, db); axis([0 1 -90 10]); grid on;
set(gca,'YTickMode','manual','YTick',[-53,-9,0])
set(gca,'YTickLabelMode','manual','YTickLabel',['53';' 9';' 0']);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.33,0.4,0.6,0.66,0.75,1]);
xlabel('frequency in pi units'); ylabel('Decibels'); title('Magnitude Response in dB');

subplot(2,2,3); plot(ww/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',[0,0.25,0.33,0.4,0.6,0.66,0.75,1])
set(gca,'YTickMode','manual','YTick',[0,1]);

subplot(2,2,4); 
sb1w = ww(1:1:ws1i)/pi; sb1e = Hr(1:1:ws1i); 
pbw = ww(wp1i:wp2i)/pi; pbe = Hr(wp1i:wp2i)-1;
sb2w = ww(ws2i:501)/pi; sb2e = Hr(ws2i:501); 
plot(sb1w,sb1e,pbw,pbe*(delta2/delta1),sb2w,sb2e);   % weighted error 
%plot(sb1w,sb1e,pbw,pbe,sb2w,sb2e);   %  error 

axis([0, 1, -deltaL, deltaL]); grid on;
xlabel('frequency in pi units'); ylabel('Hr(w)'); 
title('Weighted Error');
%title('Error Response');
set(gca,'XTickMode','manual','XTick',f)
set(gca,'YTickMode','manual','YTick',[-deltaL, 0,deltaL]);
set(gca,'XGrid','on','YGrid','on')



figure('NumberTitle', 'off', 'Name', 'Problem 7.29a Parks-McClellan Method')
set(gcf,'Color','white'); 
subplot(2,2,1); plot(w/pi, db); grid on; axis([0 2 -90 10]); 
set(gca,'YTickMode','manual','YTick',[-53,-9,0])
set(gca,'YTickLabelMode','manual','YTickLabel',['53';' 9';' 0']);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,0.6,0.75,1,1.25,1.4,1.6,1.75,2]);
xlabel('frequency in pi units'); ylabel('Decibels'); title('Magnitude Response in dB');

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');
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,0.6,0.75,1,1.25,1.4,1.6,1.75,2]);
set(gca,'YTickMode','manual','YTick',[0,1.0]);

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.29 AmpRes of h(n), Parks-McClellan Method')
set(gcf,'Color','white'); 

plot(ww/pi, Hr); grid on; %axis([0 1 -100 10]); 
xlabel('frequency in pi units'); ylabel('Hr'); title('Amplitude Response');
set(gca,'YTickMode','manual','YTick',[-delta2_db ,0,delta2_db , 1-delta1_db, 1, 1+delta1_db]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,0.6,0.75,1]);

　　运行结果：

       阶数N=22时，阻带衰减43dB，不满足要求。

       所以增大N，当N=24时，As=52dB，滤波器长度M=25