設原始信号為
采用IIR濾波器和FIR濾波器設計低通濾波器，比較兩類濾波器的濾波結果。
要求：采用MATLAB語言實作，并分析結果。
- 設計過程
首先對原始信号進行頻譜分析，确定濾波器參數。通過分析加噪信号的頻譜，噪聲信号為
，
原始信号為
。是以确定參數ωp=0.1π，ωs=0.25π；αp=0.08，αs=8；采樣間隔T=8s。
- IIR濾波器設計
按照數字濾波器技術名額（通帶邊界頻率Wp 通帶最大衰減 阻帶截止頻率 阻帶最小衰減）要求設計一個過渡模拟低通濾波器Ha(s)，再使用脈沖響應不變法或雙線性變換法将Ha(s)轉換成數字低通濾波器的系統函數H(z)。本題采用雙線性變換法，具體轉換關系為：
由參數設計巴特沃斯濾波器，根據[N,wc]=buttord(wp,ws,Rp,As,’s’)計算模拟濾波器的階數及截止頻率，由[B,A]=butter(N,wc,Rp,’s’)計算模拟濾波器系統函數的分子和分母多項式的系數向量B和A。
- FIR濾波器設計
（1）根據對阻帶衰減及過渡帶的名額要求，選擇窗函數的類型，并估計視窗的長度N。
（2）構造希望逼近的頻率響應函數Hd（e^jw）
（3）計算hd(n)
（4）加窗得到設計結果：h（n）=hd(n)*w（n）
采用[n,Wn,beta,type]=kaiserord(f,a,dev)函數來估計濾波器的階數，根據阻帶衰減及過渡帶的名額要求，選擇窗函數的類型（本題使用方法二及凱塞窗濾波器） 。這裡f對應的頻率，a=[1,0]為f指定的各個頻帶上的幅值向量，一般隻有0和1表示。Devs用于指定各個頻帶輸出濾波器的頻率響應與其期望幅值之間的最大輸出誤差或者偏差。
結果分析：
A.相同技術名額下濾波器階數
IIR濾波器階數：6；
FIR濾波器階數：14。
由此可以得到，為了達到相同的技術名額，即通帶截止頻率、阻帶起始頻率、通帶内衰減、阻帶内衰減，FIR濾波器所需的階數更高,因為FIR濾波器傳輸函數的極點固定在原點，是以隻能用較高的階數達到高的選擇性；對于同樣的濾波器設計名額，FIR濾波器所要求的階數可以比IIR濾波器高5~10倍，結果，成本較高，信号延時也較大；如果按相同的選擇性和相同的線性要求來說，則IIR濾波器就必須加全通網絡進行相位較正，同樣要大增加濾波器的節數和複雜性。
B.濾波後信号頻譜圖的分析
由圖6，FIR濾波器主要采用非遞歸結構，不論在理論上還是在實際的有限精度運算中都不存在穩定性問題，運算誤差也較小。此外，FIR濾波器可以采用快速付裡葉變換算法，在相同階數的條件下，運算速度可以快得多。可以看出，經兩種濾波器濾波後輸出的信号頻譜，FIR濾波後信号的頻譜更接近實際無噪聲信号的頻譜。
- 程式
clear;
close all;
clc
n=1:8:8000;%以8為間隔取1000個點
s(n)=cos(n/100);%噪聲信号
x(n)=cos(n/100)+0.01*sin(n);%原始信号
n1=1:512;%對采樣後序列進行512點傅裡葉變換
fk=1:8000/512:8000;%橫坐标轉化為頻率
%無噪信号和原始信号
figure(1)
subplot(1,2,1)
plot(n,s(n))
xlabel('n');
ylabel('s(n)');
title('無噪信号');
hold on
subplot(1,2,2)
plot(n,x(n))
xlabel('n');
ylabel('x(n)');
title('原始信号');
fw=fft(s(n),512);
fw1=fft(x(n),512);
mag=abs(fw);
mag1=abs(fw1);
%原始信号頻譜分析
figure(2)
subplot(1,2,1)
plot(fk/8000,20*log10(mag),'b')
xlabel('w');
ylabel('幅度');
title('無噪信号頻譜');
hold on
subplot(1,2,2)
plot(fk/8000,20*log10(mag1),'r')
xlabel('w');
ylabel('幅度');
title('原始信号頻譜');
%IIR濾波器設計
wp=0.1*pi;%通帶截止頻率
ws=0.25*pi; %阻帶起始頻率
Ap=1; %通帶内衰減
As=8; %阻帶衰減
T=1;
WP=(2/T)*tan(wp/2);
WS=(2/T)*tan(ws/2);
[N,Wc]= buttord(WP,WS,Ap,As,'s');%傳回濾波器階數和3dB截止頻率
[bt,at]=butter(N,Wc,'s');%計算低通濾波器系統函數H(s)
[bd,ad ]=bilinear(bt,at,1/T); %用雙線性變換法設計數字濾波器，數字濾波器的頻率響應
[hd,wd]=freqz(bd,ad,512); %畫圖
fw2=fft(hd,512);
mag2=abs(fw2);
figure
plot(fk/8000,20*log10(mag2),'r')
hold on
plot(fk/8000,20*log10(mag1),'b')
title('原始信号頻譜和IIR濾波器頻譜對比圖');
legend('IIR濾波器頻譜','原始信号頻譜')
y=filter(bd,ad,x(1));
%%FIR濾波器的設計
fcuts = [0.1 0.25]; %歸一化頻率omega/pi，這裡指通帶截止頻率、阻帶起始頻率
mags = [1 0];
devs = [10^(-Ap/20) (10^(As/20)-1)/(10^(As/20+1))];
[nn,Wn,beta,ftype] = kaiserord(fcuts,mags,devs); %計算出凱塞窗N，beta的值
hh=kaiser(nn+1,beta);
fw2=fft(hh,512);
mag22=abs(fw2);
figure
plot(fk/8000,20*log10(mag1),'r')
hold on
plot(fk/8000,20*log10(mag22),'b')
xlabel('w');
ylabel('幅度');
title('原始信号頻譜和FIR濾波器頻譜對比圖');
legend('原始信号頻譜','FIR濾波器頻譜');
%去噪效果對比
figure
y1=filter(hh,5,x(n));
plot(y1,'g')
hold on
plot(y,'r')
hold on
plot(s(n),'b')
title('原始信号頻譜和濾波器頻譜對比圖');
legend('FIR','IIR','無噪信号')
figure
fwww=fft(y,512);
magnewiir=abs(fwww);
plot(fk/4000,20*log(magnewiir),'r')
hold on
fwwww=fft(y1,512);
magnewfir=abs(fwwww);
plot(fk/4000,magnewfir,'b')
hold on
plot(fk/4000,20*log(mag1),'k')
legend('iir濾波後頻譜','fir濾波後頻譜','原始有噪聲信号頻譜');