《DSP using MATLAB》示例 Example 9.11

Posted 2020-09-20 沧海一粟

代码：

%% ------------------------------------------------------------------------
%%            Output Info about this m-file
fprintf(\'\\n***********************************************************\\n\');
fprintf(\'        <DSP using MATLAB> Exameple 9.11 \\n\\n\');

time_stamp = datestr(now, 31);
[wkd1, wkd2] = weekday(today, \'long\');
fprintf(\'      Now is %20s, and it is %7s  \\n\\n\', time_stamp, wkd2);
%% ------------------------------------------------------------------------

% Filter Design:
D = 2; Rp = 0.1; As = 30; wp = pi/D; ws = wp+0.1*pi; 
[delta1, delta2] = db2delta(Rp, As); 
%weights = [delta2/delta1, 1];
[N, F, A, weights] = firpmord([wp, ws]/pi, [1, 0], [delta1, delta2], 2); 
h = firpm(N, F, A, weights); n = [0:length(h)-1];
[Hr, w, a, L] = Hr_Type1(h); Hr_min = min(Hr); w_min = find(Hr == Hr_min);
H = abs(freqz(h, 1, w)); Hdb = 20*log10(H/max(H)); min_attn = Hdb(w_min);

%% -----------------------------------------------------------------
%%                             Plot
%% -----------------------------------------------------------------  

% Input signal
Hf1 = figure(\'units\', \'inches\', \'position\', [1, 1, 8, 6], ...
	\'paperunits\', \'inches\', \'paperposition\', [0, 0, 6, 4], ...
	\'NumberTitle\', \'off\', \'Name\', \'Exameple 9.11\');
set(gcf,\'Color\',\'white\'); 

TF = 10;

subplot(2, 2, 1); 
Hs1 = stem(n, h, \'filled\'); set(Hs1, \'markersize\', 2, \'color\', \'g\');
axis([-1, length(n), -0.15, 0.6]); grid on;
xlabel(\'n\', \'vertical\', \'bottom\'); ylabel(\'Amplitude\', \'vertical\', \'cap\'); 
title(\'Inpulse Response \', \'fontsize\', TF, \'vertical\', \'baseline\');
set(gca, \'xtick\', [n(1), n(end)]); 
set(gca, \'ytick\', [0, 0.5]);


subplot(2, 2, 3); 
plot(w/pi, Hr, \'m\', \'linewidth\', 1.0); axis([0, 1, -0.1, 1.1]); grid on;
xlabel(\'Frequency in \\pi units\', \'vertical\', \'middle\'); ylabel(\'Amplitude\', \'vertical\', \'cap\'); 
title(\'Amplitude Response\', \'fontsize\', TF, \'vertical\', \'baseline\');
set(gca, \'xtick\', [0, wp/pi, ws/pi, 1]); 
set(gca, \'ytick\', [0, 1]);


subplot(2, 2, 2);
plot(w/pi, Hdb, \'m\', \'linewidth\', 1.0); axis([0, 1, -50, 10]); grid on;
xlabel(\'Frequency in \\pi units\', \'vertical\', \'middle\'); ylabel(\'Decibels\', \'vertical\', \'cap\'); 
title(\'Log-magnitude Response \', \'fontsize\', TF, \'vertical\', \'baseline\');
set(gca, \'xtick\', [0, wp/pi, ws/pi, 1]); 
set(gca, \'ytick\', [-50, round(min_attn), 0]);

subplot(2, 2, 4); 
lw = length(w)-1; PB = [0:floor(wp/pi*lw)]; HrPB = Hr(PB+1)-1;
SB = [ceil(ws/pi*lw):lw]; HrSB = Hr(SB+1);
[AX, H1, H2] = plotyy(PB/lw, HrPB, SB/lw, HrSB);
delta1 = round(delta1*1000)/1000; delta2 = round(delta2*100)/100;
set(AX(1), \'xtick\', [0, wp/pi, ws/pi, 1], \'ytick\', [-delta1, 0, delta1], \'Ycolor\', \'g\');
set(AX(2), \'xtick\', [0, wp/pi, ws/pi, 1], \'ytick\', [-delta2, 0, delta2], \'Ycolor\', \'r\');

set(H1, \'color\', \'g\', \'linewidth\', 1);
set(H2, \'color\', \'r\', \'linewidth\', 1);
title(\'Unweighted Ripples\', \'fontsize\', TF, \'vertical\', \'baseline\');

xlabel(\'Frequency in \\pi units\', \'vertical\', \'middle\'); ylabel(\'Amplitude\', \'vertical\', \'cap\'); 

　　运行结果：

    滤波器在通带[0，π/2]上无畸变通过信号。但是，因为过渡带较宽，所以过渡带中的信号有可能假频混入到感兴趣的频带
中。30dB的衰减可能使得在减采样过程后阻带的部分谱折叠到通带中。因此，我们需要一个更好的方法。