測試代碼：

%%****bin_bp.m*******************%%
%% 使用雙線性變換法設計帶通濾波器
%% 2018年6月13日 16:30:34
%% author:Alimy

close all;
clear;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%代碼正文
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%給定數字濾波器指標
f_sl    =   150     ;   %阻帶下限頻率(Hz)
f_1     =   200     ;   %通帶下限頻率(Hz)
f_3     =   500     ;   %通帶上限頻率(Hz)
f_sh    =   600     ;   %阻帶上限頻率(Hz)
R_p     =   0.5     ;   %通帶允許的最大衰減
R_s     =   40      ;   %阻帶允許的最小衰減
f_s     =   2000    ;   %采樣頻率
T_s     =   1 / f_s ;   %采樣間隔
%1.將數字帶通濾波器的頻率參數變換為歸一化的數字角頻率參數
omega_sl    =   2 * pi * f_sl   / f_s;    %阻帶下限頻率
omega_1     =   2 * pi * f_1    / f_s;    %通帶下限頻率
omgea_3     =   2 * pi * f_3    / f_s;    %通帶上限頻率
omega_sh    =   2 * pi * f_sh   / f_s;    %阻帶上限頻率
%2.預畸變處理,將歸一化數字角頻率參數變換成模擬帶通濾波器的角頻率參數
C = 2*f_s ;
Omega_sl    =   C * tan( omega_sl  / 2 );
Omega_1     =   C * tan( omega_1  / 2 );
Omega_3     =   C * tan( omgea_3  / 2 );
Omega_sh    =   C * tan( omega_sh  / 2 );
%3.對模擬帶通濾波器的角頻率參數做歸一化處理
Omega_BW    =   Omega_3 - Omega_1;
eta_sl      =   Omega_sl / Omega_BW;
eta_1       =   Omega_1  / Omega_BW;   
eta_3       =   Omega_3  / Omega_BW;
eta_sh      =   Omega_sh / Omega_BW;
%4.設計歸一化模擬濾波器,得到歸一化模擬帶通濾波器的轉移函數
Omega_p = [ Omega_1  , Omega_3 ];
Omega_s = [ Omega_sl , Omega_sh ];
[ N , Wn ] = buttord( Omega_p , Omega_s , R_p , R_s , ‘s‘ ); %選擇模擬巴特沃斯濾波器的最小階數
[ Bs, As ] = butter(N,Wn,‘s‘);
%5.利用模擬帶通濾波器的轉移函數確定IIR數字濾波器的轉移函數 
[ bz , az ] = bilinear(Bs,As,f_s);

figure(1);
freqz(bz,az);
title(‘帶通濾波器幅度譜和相位譜特性‘);
%濾波效果測試 
N = 1000;
t = [ 0 : N - 1 ] * T_s ;
f1 = 50;
f2 = 100;
f3 = 300;
x1 = 2*1*sin( 2 * pi * f1 * t );
x2 = 2*2*sin( 2 * pi * f2 * t );
x3 = 2*1*sin( 2 * pi * f3 * t );
x  = x1 + x2 + x3;

fft_x = fft( x );
X_mag = fftshift( abs   ( fft_x ) ) / N ;
X_ang = fftshift( angle ( fft_x ) );
delta_f = f_s/N;
f = ( -N / 2 : N / 2 - 1 )*delta_f;

figure(2);
subplot(3,1,1);
plot(t,x);
title(‘原信號時域波形‘);
xlabel(‘t(s)‘);
subplot(3,1,2);
plot( f , X_mag );
title(‘原信號幅度譜‘);
xlabel(‘f(hz)‘);
subplot(3,1,3);
plot( f , X_ang );
title(‘原信號相位譜‘);
xlabel(‘f(hz)‘);


%濾波
lp_x = filter( bz , az , x );
lp_fft_x = fft( lp_x );
lp_X_mag = fftshift( abs   ( lp_fft_x ) ) / N ;
lp_X_ang = fftshift( angle ( lp_fft_x ) );
figure(3);
subplot(3,1,1);
plot(t,lp_x);
title(‘濾波後信號時域波形‘);
xlabel(‘t(s)‘);
subplot(3,1,2);
plot( f , lp_X_mag );
title(‘濾波後信號幅度譜‘);
xlabel(‘f(hz)‘);
subplot(3,1,3);
plot( f , lp_X_ang );
title(‘濾波後信號相位譜‘);
xlabel(‘f(hz)‘);

　效果：

濾波器特性：

待濾波的信號：

濾波後的信號:

[Matlab]雙線性變換法設計數字帶通濾波器