上代码：
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('\n***********************************************************\n');
fprintf(' <DSP using MATLAB> Problem 8.36.1 \n\n');
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('\n*******Digital lowpass, Coefficients of DIRECT-form***********\n');
[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('\nActual Passband Ripple is %.4f dB.\n', Rp_lp);
As_lp = -round(max(dblp(ceil(wslp/delta_w)+1:1:501))); % Min Stopband attenuation
fprintf('\nMin Stopband attenuation is %.4f dB.\n\n', 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 = '\nPlease input the type of digital lp filter: \n\n 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('\nActual Passband Ripple is %.4f dB.\n', Rp_lp);
As_lp = -round(max(dblp(ceil(wslp/delta_w)+1 :1 : 501))); % Min Stopband attenuation
fprintf('\nMin Stopband attenuation is %.4f dB.\n\n', 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
可见自带函数比个人所写的效果强！