%巴特沃斯數字濾波器設計
clc;
close?all;
clear?all;

%%?巴特沃斯濾波器
%原始信號

fs=30;??????????????%采樣頻率
N=300;??????????????%N/fs?秒數據
n=0:N-1;????????????
t=n/fs;?????????????%時間
%信號形式?
fl=0.4;???????????%頻率1（低頻）
f2=2;?????????????%頻率2
f3=5;?????????????%頻率3（高頻）
s=cos（2*pi*fl*t）+cos（2*pi*f2*t）+cos（2*pi*f3*t）;????%s是0.4Hz、2Hz和5Hz信號疊加，低通截止頻率是1Hz
subplot（121）;plot（ts）;
title（‘輸入信號‘）;xlabel（‘t/s‘）;ylabel（‘幅度‘）;
sfft=fft（s）;
subplot（122）;
plot（（1:length（sfft）/2）*fs/length（sfft）2*abs（sfft（1:length（sfft）/2））/length（sfft））;
title（‘信號頻譜‘）;xlabel（‘頻率/Hz‘）;ylabel（‘幅度‘）;

%%?設計低通濾波器，截止頻率為1

Wp_l=1/（fs/2）;?%Wp：通帶截止頻率為1Hz
Ws_l=2/（fs/2）;?%Ws：阻帶截止頻率為2Hz
%估算得到Butterworth低通濾波器的最小階數N和3dB截止頻率Wn
Rp_l=1;%通帶最大衰減Rp=1，通帶紋波小于1db
Rs_l=50;%阻帶最小衰減Rs=50（dB），阻帶衰減大于50db
[n_lWn_l]=buttord（Wp_lWs_lRp_lRs_l）;%

%設計Butterworth低通濾波器
[a_lb_l]=butter（n_lWn_l）;%默認為低通和帶通濾波器
[h_lf_l]=freqz（a_lb_l‘whole‘fs）;????????%求數字低通濾波器的頻率響應
f_l=（0:length（f_l）-1*fs/length（f_l））;?????%進行對應的頻率轉換
figure（2）;
plot（f_l（1:length（f_l）/2）abs（h_l（1:length（f_l）/2）））;???????%繪制幅頻響應圖
title（‘巴特沃斯低通濾波器‘）;xlabel（‘頻率/Hz‘）;ylabel（‘幅度‘）;
grid;
sF_l=filter（a_lb_ls）;???????????????????%疊加函數s經過低通濾波器以后的新函數
figure（3）;
subplot（121）;
plot（tsF_l）;?????????????????????????%繪制疊加函數s經過低通后時域圖形
title（‘輸出信號‘）;xlabel（‘t/s‘）;ylabel（‘幅度‘）;
SF_l=fft（sF_l）;
subplot（122）;
plot（（1:length（SF_l）/2）*fs/length（SF_l）2*abs（SF_l（1:length（SF_l）/2））/length（SF_l））;
title（‘低通濾波后頻譜‘）;xlabel（‘頻率/Hz‘）;ylabel（‘幅度‘）;

%%?高通
%?wn=1/（fs/2）;

%?Wp_h=4/（fs/2）;%0.6*pi;%通帶拐角頻率Wp屬于[0，1]，Wp=1時對應于歸一化角頻率pi
%?Ws_h=3/（fs/2）;%0.4*pi;%阻帶起始頻率Wp屬于[0，1]
%?Rp_h=1;%通帶內最大衰減，通帶內波動（db）
%?Rs_h=15;%阻帶內最小衰減（db）
%?[N_hwn_h]=buttord（Wp_h/piWs_h/piRp_hRs_h）;?%計算巴特沃斯濾波器階次和截止頻率Wn
wn_h=3/（fs/2）;%截止頻率3HZ
N_h=10;
[a_hb_h]=butter（N_hwn_h‘high‘）;??????????%頻率變換法設計巴特沃斯高通濾波器
[Hw_h]=freqz（a_hb_hfs）;???????????????????????%500點的復頻響應
mag=abs（H）;?????????????????????????????????%絕對幅值響應
%?db=20*log10（mag/max（mag））;??????????????????%相對幅值響應