實驗目的和要求
1. 設計IIR巴特沃斯低通濾波器
2. 將濾波器用於加噪聲訊號的處理
3. 改變引數指標比較不同
4.
實驗結果與分析
（1）濾波器頻率特性如圖

輸入訊號特性如圖:

經過濾波器之後訊號特性如圖:

（2）改變指標後
Wp=0.3*pi, Ws=0.6*pi, Rp=0.5, As=40

%課本P395
 

%1、由模擬濾波器指標（Omegap、Omegas、Rp、As）計算巴特沃斯濾波器指標（階數 N 和 3dB 截止頻率Omegac）
%buttord可用於計算低通、高通、帶通、帶阻巴特沃斯濾波器指標
%Rp，As分別為通帶最大衰減和組帶最小衰減（dB），'s'表示模擬濾波器
%根據指標，求得符合指標的、濾波器應有的階數和截止頻率
clc;
Td=1;Fs=1/Td;
Wp=0.2*pi;Ws=0.3*pi;
Rp=1;
As=15;
Omegap=2/Td*tan(Wp/2);
Omegas=2/Td*tan(Ws/2);
[N,Omegac]=buttord(Omegap,Omegas,Rp,As,'s'
 
);
%2、計算連續系統係數a,b0
%歸一化（模板巴特沃斯的截止頻率Omegac=1）的模擬濾波器，零點、極點、增益
%級聯式，z0,p0,k0分別為零點、極點、增益
%去歸一化（求出真正濾波器的Omegac）,得到級聯式
%級聯式轉換為一般式（由根求係數）：為由一般式為跳板，計算並聯式
[z0,p0,k0]=buttap(N);
[Bap,Aap]=zp2tf(z0,p0,k0); 
[b,a]=lp2lp(Bap,Aap,Omegac);
[bz,az]=bilinear(b,a,Fs); 
[H,W]=freqs(bz,az);
figure;
plot(W/pi,20*log10(abs(H))); 
title('巴特沃斯濾波器'
 
)
xlabel('頻率');
ylabel('增益(dB)');

f1=10;%第一個點頻訊號分量頻率
f2=30;%第二個點頻訊號分量頻率
f3=45;%第三個點頻訊號分量頻率
fs=100;%取樣率
T=2;%間隔長度
n=round(T*fs);%取樣點個數
t=linspace(0,T,n);
y=cos(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t)+randn(size(t));
figure;
subplot(2,1,1);
plot(t,y);
title('輸入訊號時域影象');
xlabel('t/s');
ylabel('V');
fft_y=fftshift(fft(y));
f=linspace(-fs/2,fs/2,n);
subplot(2,1,2);
plot(f,abs(fft_y));
title('輸入訊號頻域影象');
xlabel('f/Hz');
ylabel('V');
axis([ 0 50 0 100]);

figure;
final=filter(bz, az, y); 
subplot(2,1,1);
plot(t,final);  
title( '濾波後訊號波形時域波形' );
xlabel('t/s');
ylabel('V');

finall=fftshift(fft(final));
subplot(2,1,2);
plot(f,abs(finall));
title( '濾波後訊號頻域波形' );
axis([ 0 50 0 100]);

注：課本指程佩青，數字訊號處理教程，清華大學出版社。