《DSP using MATLAB》Problem 8.39

代码：

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

');

banner();
%% ------------------------------------------------------------------------

% Digital Filter Specifications:   Chebyshev-1 bandstop
ws = [0.35*pi 0.65*pi];             % digital stopband freq in rad
wp = [0.25*pi 0.75*pi];             % digital passband freq in rad

delta1 = 0.05;                      % passband tolerance, absolute specs 
delta2 = 0.01;                      % stopband tolerance, absolute specs

Rp = -20 * log10( (1-delta1)/(1+delta1))    % passband ripple in dB
As = -20 * log10( delta2 / (1+delta1))      % stopband attenuation in dB                       

Ripple = 10 ^ (-Rp/20)           % passband ripple in absolute
Attn = 10 ^ (-As/20)             % stopband attenuation in absolute

% Calculation of Chebyshev-1 filter parameters:
[N, wn] = cheb1ord(wp/pi, ws/pi, Rp, As);

fprintf('
  ********* Chebyshev-1 Digital Bandstop Filter Order is = %3.0f 
', 2*N)

% Digital Chebyshev-1 bandstop Filter Design:
[bbs, abs] = cheby1(N, Rp, wn, 'stop');

[C, B, A] = dir2cas(bbs, abs)

% Calculation of Frequency Response:
[dbbs, magbs, phabs, grdbs, wwbs] = freqz_m(bbs, abs);

% ---------------------------------------------------------------
%    find Actual Passband Ripple and Min Stopband attenuation
% ---------------------------------------------------------------
delta_w = 2*pi/1000;
Rp_bs = -(min(dbbs(1:1:ceil(wp(1)/delta_w+1))));                           % Actual Passband Ripple

fprintf('
Actual Passband Ripple is %.4f dB.
', Rp_bs);

As_bs = -round(max(dbbs(ceil(ws(1)/delta_w)+1:1:ceil(ws(2)/delta_w)+1)));   % Min Stopband attenuation
fprintf('
Min Stopband attenuation is %.4f dB.

', As_bs);


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

figure('NumberTitle', 'off', 'Name', 'Problem 8.39.2 Chebyshev-1 bs by cheby1 function')
set(gcf,'Color','white'); 
M = 1;                          % Omega max

subplot(2,2,1); plot(wwbs/pi, magbs); axis([0, M, 0, 1.2]); grid on;
xlabel('Digital frequency in pi units'); ylabel('|H|'); title('Magnitude Response');
set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.25, 0.35, 0.65, 0.75, M]);
set(gca, 'YTickMode', 'manual', 'YTick', [0, 0.01, 0.9048, 1]);

subplot(2,2,2); plot(wwbs/pi, dbbs); axis([0, M, -120, 2]); grid on;
xlabel('Digital frequency in pi units'); ylabel('Decibels'); title('Magnitude in dB');
set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.25, 0.35, 0.65, 0.75, M]);
set(gca, 'YTickMode', 'manual', 'YTick', [-80, -40, -1, 0]);
set(gca,'YTickLabelMode','manual','YTickLabel',['80'; '40';'1 ';' 0']);


subplot(2,2,3); plot(wwbs/pi, phabs/pi); axis([0, M, -1.1, 1.1]); grid on;
xlabel('Digital frequency in pi nuits'); ylabel('radians in pi units'); title('Phase Response');
set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.25, 0.35, 0.65, 0.75, M]);
set(gca, 'YTickMode', 'manual', 'YTick', [-1:0.5:1]);

subplot(2,2,4); plot(wwbs/pi, grdbs); axis([0, M, 0, 50]); grid on;
xlabel('Digital frequency in pi units'); ylabel('Samples'); title('Group Delay');
set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.25, 0.35, 0.65, 0.75, M]);
set(gca, 'YTickMode', 'manual', 'YTick', [0:20:50]);

　　运行结果：

        通带、阻带衰减指标，dB单位和绝对值单位，

        使用cheby1函数（MATLAB工具箱函数）得到Chebyshev-1型数字带阻滤波器，系统函数串联形式的系数

       最小阻带衰减为44dB，满足40.4238dB的设计要求。滤波器的幅度谱、相位谱和群延迟响应如下

        类似的，使用butter、cheby2和ellip函数（都是MATLAB工具箱函数）也可以设计带阻滤波器，这里不放图了。