《DSP using MATLAB》Problem 8.36

上代码：

function [wpLP, wsLP, alpha] = lp2lpfre(wplp, wslp)
% Band-edge frequency conversion from lowpass to lowpass digital filter
% -------------------------------------------------------------------------
% [wpLP, wsLP, alpha] = lp2lpfre(wplp, wslp) 
%   wpLP = passband edge for the lowpass digital prototype 
%   wsLP = stopband edge for the lowpass digital prototype 
%  alpha = lowpass to lowpass transformation parameter
%   wplp = passband edge frequency for the given lowpass
%   wslp = stopband edge frequency for the given lowpass
%
%
if wplp <= 0 
	error('Passband edge must be larger than 0.')
end

if wslp <= wplp
	error('Stopband edge must be larger then Passband edge')
end

% Determine the digital lowpass cutoff frequencies:
wpLP = 0.2*pi;
alpha = sin((wpLP - wplp)/2)/sin((wpLP + wplp)/2);
wsLP = -angle((exp(-j*wslp)-alpha)/(1-alpha*exp(-j*wslp)));

　　

function [b, a] = dlpfd_bl(type, wp, ws, Rp, As)
% IIR Lowpass Filter Design using bilinear transformation
% -----------------------------------------------------------------------
% [b, a] = dlpfd_bl(type, wp, ws, Rp, As);
%   type = 'butter' or 'cheby1' or 'cheby2' or 'ellip'   
%      b = numerator polynomial coefficients of lowpass filter , Direct form
%      a = denominator polynomial coefficients of lowpass filter, Direct form
%     wp = Passband edge frequency in radians;
%     ws = Stopband edge frequency in radians (wp < ws);
%     Rp = Passband ripple in +dB; Rp > 0
%     As = Stopband attenuation in +dB; As > 0
%
%
%prompt = 'Please input the type of digital lp filter:  butter or cheby1 or cheby2 or ellip [butter] ';
%type = input(prompt , 's');

if isempty(type)
    str = 'butter';
end

switch type
    case 'butter'
    	[b , a] = buttlpf(wp, ws, Rp, As);
    case 'cheby1'
    	[b , a] = cheb1lpf(wp, ws, Rp, As);
    case 'cheby2'
    	[b , a] = cheb2lpf(wp, ws, Rp, As);
    case 'ellip'
    	[b , a] = eliplpf(wp, ws, Rp, As);
    otherwise
        disp('Oh, input may be error!');
end

　　第1小题

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

');

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

% Digital lowpass Filter Specifications:
wplp = 0.45*pi;                 % digital passband freq in rad
wslp = 0.50*pi;                 % digital stopband freq in rad
Rp = 0.5;                       % passband ripple in dB
As = 60;                        % stopband attenuation in dB

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

fprintf('
*******Digital lowpass, Coefficients of DIRECT-form***********
');
[blp, alp] = buttlpf(wplp, wslp, Rp, As)
%[blp, alp] = cheb1lpf(wphp, wshp, Rp, As)
%[blp, alp] = cheb2lpf(wphp, wshp, Rp, As)
%[blp, alp] = eliplpf(wphp, wshp, Rp, As)
[C, B, A] = dir2cas(blp, alp)

% Calculation of Frequency Response:
[dblp, maglp, phalp, grdlp, wwlp] = freqz_m(blp, alp);
%[dbhp, maghp, phahp, grdhp, wwhp] = freqz_m(bhp, ahp);

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

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

As_lp = -round(max(dblp(ceil(wslp/delta_w)+1:1:501)));                    % Min Stopband attenuation
fprintf('
Min Stopband attenuation is %.4f dB.

', As_lp);

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

figure('NumberTitle', 'off', 'Name', 'Problem 8.36.1 Butterworth lowpass by buttlpf function')
set(gcf,'Color','white'); 
M = 2;                          % Omega max

subplot(2,2,1); plot(wwlp/pi, maglp); axis([0, M, 0, 1.2]); grid on;
xlabel('Digital frequency in pi units'); ylabel('|H|'); title('Lowpass Filter Magnitude Response');
set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.45, 0.5, M]);
set(gca, 'YTickMode', 'manual', 'YTick', [0, 0.9441, 1]);

subplot(2,2,2); plot(wwlp/pi, dblp); axis([0, M, -150, 1]); grid on;
xlabel('Digital frequency in pi units'); ylabel('Decibels'); title('Lowpass Filter Magnitude in dB');
set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.45, 0.5, M]);
set(gca, 'YTickMode', 'manual', 'YTick', [ -111, -85, -1, 0]);
%set(gca,'YTickLabelMode','manual','YTickLabel',['111'; '85'; '1 ';' 0']);


subplot(2,2,3); plot(wwlp/pi, phalp/pi); axis([0, M, -1.1, 1.1]); grid on;
xlabel('Digital frequency in pi nuits'); ylabel('radians in pi units'); title('Lowpass Filter Phase Response');
set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.45, 0.5, M]);
set(gca, 'YTickMode', 'manual', 'YTick', [-1:1:1]);

subplot(2,2,4); plot(wwlp/pi, grdlp); axis([0, M, 0, 20]); grid on;
xlabel('Digital frequency in pi units'); ylabel('Samples'); title('Lowpass Filter Group Delay');
set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.45, 0.5, M]);
set(gca, 'YTickMode', 'manual', 'YTick', [0:10:20]);



% -----------------------------------------------------
%              method 2
% -----------------------------------------------------
% Digital lowpass Filter Specifications:
[wpLP, wsLP, alpha] = lp2lpfre(wplp, wslp);

prompt = '
Please input the type of digital lp filter: 

  butter or cheby1 or cheby2 or ellip [butter]: ';
type = input(prompt , 's');

[blp, alp] = dlpfd_bl(type, wplp, wslp, Rp, As);

[C, B, A] = dir2cas(blp, alp);



% -----------------------------------------------------
%              method 3  butter function
% -----------------------------------------------------
% Calculation of Butterworth lp filter parameters:
[N, wn] = buttord(wplp/pi, wslp/pi, Rp, As)

% Digital Butterworth lowpass Filter Design:
[blp, alp] = butter(N, wn, 'low')

[C, B, A] = dir2cas(blp, alp)

% Calculation of Frequency Response:
[dblp, maglp, phalp, grdlp, wwlp] = freqz_m(blp, alp);
%[dbhp, maghp, phahp, grdhp, wwhp] = freqz_m(bhp, ahp);

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

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

As_lp = -round(max(dblp(ceil(wslp/delta_w)+1 :1 : 501)));                    % Min Stopband attenuation
fprintf('
Min Stopband attenuation is %.4f dB.

', As_lp);

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

figure('NumberTitle', 'off', 'Name', 'Problem 8.36.1 Butterworth lowpass by butter function')
set(gcf,'Color','white'); 
M = 1;                          % Omega max

subplot(2,2,1); plot(wwlp/pi, maglp); axis([0, M, 0, 1.2]); grid on;
xlabel('Digital frequency in pi units'); ylabel('|H|'); title('Lowpass Filter Magnitude Response');
set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.45, 0.5, M]);
set(gca, 'YTickMode', 'manual', 'YTick', [0, 0.9441, 1]);

subplot(2,2,2); plot(wwlp/pi, dblp); axis([0, M, -100, 2]); grid on;
xlabel('Digital frequency in pi units'); ylabel('Decibels'); title('Lowpass Filter Magnitude in dB');
set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.45, 0.5, M]);
set(gca, 'YTickMode', 'manual', 'YTick', [-70, -60, -1, 0]);
set(gca,'YTickLabelMode','manual','YTickLabel',['70'; '60';'1 ';' 0']);


subplot(2,2,3); plot(wwlp/pi, phalp/pi); axis([0, M, -1.1, 1.1]); grid on;
xlabel('Digital frequency in pi nuits'); ylabel('radians in pi units'); title('Lowpass Filter Phase Response');
set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.45, 0.5, M]);
set(gca, 'YTickMode', 'manual', 'YTick', [-1:1:1]);

subplot(2,2,4); plot(wwlp/pi, grdlp); axis([0, M, 0, 90]); grid on;
xlabel('Digital frequency in pi units'); ylabel('Samples'); title('Lowpass Filter Group Delay');
set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.45, 0.5, M]);
set(gca, 'YTickMode', 'manual', 'YTick', [0:10:90]);

　　运行结果：

        绝对指标

          数字低通，频带边界截止频率

        采用buttlpf函数，数字低通butterworth滤波器阶数51，系统函数直接形式系数

        转换成串联形式，系数

        采用butter（MATLAB自带函数），计算数字低通滤波器，阶数51

        可见自带函数比个人所写的效果强！

         第3小题，Chebyshev-2型

        

        采用cheb2lpf函数，得到的Chebyshev-2型数字低通滤波器，幅度谱、相位谱和群延迟响应

        采用cheby2（MATLAB自带函数），计算得到数字低通滤波器，系统函数串联形式系数

        Chebyshev-1型和Elliptic型数字低通，这里不放图了。