《DSP using MATLAB》示例Example7.15
代码:
%T1 = 0.5 M = 40; alpha = (M-1)/2; l = 0:M-1; wl = (2*pi/M)*l; Hrs = [ones(1, 5), 0.5, zeros(1, 29), 0.5, ones(1, 4)]; % Ideal Amp Res sampled Hdr = [1, 1, 0, 0]; wdl = [0, 0.25, 0.25, 1]; % Ideal Amp Res for plotting k1 = 0:floor((M-1)/2); k2 = floor((M-1)/2)+1:M-1; %% -------------------------- %% Type-2 LPF %% -------------------------- 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 a figure('NumberTitle', 'off', 'Name', 'Exameple 7.15a') 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('Frequency Samples: M=40, T1=0.5'); grid on; subplot(2,2,2); stem(l, h); axis([-1, M, -0.1, 0.3]); 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, -60, 10]); grid on; xlabel('frequency in \pi units'); ylabel('Decibels'); title('Magnitude Response');
运行结果:
最小阻带衰减变成30dB,这比7.14例子要好些。但仍不符要求(50dB)。最佳的T1值需要手动修改来得到,
当然有效的方法是通过线性编程技巧来获得T1,本书中没有提及。最优解T1=0.39,如下:
%T1 = 0.39 M = 40; alpha = (M-1)/2; l = 0:M-1; wl = (2*pi/M)*l; Hrs = [ones(1, 5), 0.39, zeros(1, 29), 0.39, ones(1, 4)]; % Ideal Amp Res sampled Hdr = [1, 1, 0, 0]; wdl = [0, 0.25, 0.25, 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 b figure('NumberTitle', 'off', 'Name', 'Exameple 7.15b') 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]); grid on; xlabel('frequency in \pi nuits'); ylabel('Hr(k)'); title('Frequency Samples: M=40, T1=0.39'); subplot(2,2,2); stem(l, h); axis([-1, M, -0.1, 0.3]); 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, -60, 10]); grid on; xlabel('frequency in \pi units'); ylabel('Decibels'); title('Magnitude Response');
运行结果:
最优阻带衰减现在是43dB,很明显如果想进一步增加阻带衰减As,那么就需要变动过渡带中1个以上的样点。
牢记:
1、如果你决定做某事,那就动手去做;不要受任何人、任何事的干扰。2、这个世界并不完美,但依然值得我们去为之奋斗。