《DSP using MATLAB》示例 Example 9.11

代码：

%% ------------------------------------------------------------------------

%%            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的衰减可能使得在减采样过程后阻带的部分谱折叠到通带中。因此，我们需要一个更好的方法。