《DSP using MATLAB》示例Example7.17

代码：

M = 40; alpha = (M-1)/2; l = 0:M-1; wl = (2*pi/M)*l; 
T1 = 0.109021; T2 = 0.59417456;
Hrs = [zeros(1, 5),T1,T2, ones(1, 7),T2,T1, zeros(1, 9),T1,T2, ones(1,7), T2,T1,zeros(1,4)];    % Ideal Amp Res sampled
Hdr = [0, 0, 1, 1, 0, 0]; wdl = [0, 0.2, 0.35, 0.65, 0.8, 1];                                  % Ideal Amp Res for plotting
k1 = 0:floor((M-1)/2); k2 = floor((M-1)/2)+1:M-1;

angH = [-alpha*(2*pi)/M*k1, alpha*(2*pi)/M*(M-k2)]; 
H = Hrs.*exp(j*angH); h = real(ifft(H, M));

[db, mag, pha, grd, w] = freqz_m(h, 1);
[Hr, ww, a, L] = Hr_Type2(h);

%Plot 

figure('NumberTitle', 'off', 'Name', 'Exameple 7.17a')
set(gcf,'Color','white'); 

subplot(2,2,1); plot(wl(1:21)/pi, Hrs(1:21), 'o', wdl, Hdr); axis([0, 1, -0.1, 1.1]);
xlabel('frequency in pi nuits'); ylabel('Hr(k)'); title('Bandpass: M=40, T1=0.5941, T2=0.109');
grid on;

subplot(2,2,2); stem(l, h); axis([-1, M, -0.4, 0.4]); grid on;
xlabel('n'); ylabel('h(n)'); title('Impulse Response');

subplot(2,2,3); plot(ww/pi, Hr, wl(1:21)/pi, Hrs(1:21), 'o'); axis([0, 1, -0.2, 1.2]); grid on;
xlabel('frequency in pi units'); ylabel('Hr(w)'); title('Amplitude Response');

subplot(2,2,4); plot(w/pi, db); axis([0, 1, -100, 10]); grid on;
xlabel('frequency in pi units'); ylabel('Decibels'); title('Magnitude Response');

　　运行结果：

        从幅度谱上看出，最小阻带衰减As大概60dB，满足设计要求。