~~ 如果有什么問題可以在我的個人博客留言 ，我會及時回復。歡迎來訪交流 ~~
IIR數字濾波器設計（數字信號處理）
一、實驗目的
1.熟悉雙線性變換法設計IIR數字濾波器的原理與方法。
2.掌握IIR數字濾波器的MATLAB實現方法，會調用ellipord()和ellip()
函數設計各種濾波器。
3.觀察分析濾波器輸入輸出數據波形，理解數字濾波的概念。
二、實驗原理及步驟
（一）實驗原理-雙線性變換法
數字濾波器是對數字信號實現濾波的線性時不變系統。數字濾波實質上是一種運算過程，實現對信號的運算處理。輸入數字信號（數字序列）通過特定的運算轉變為輸出的數字序列，因此，數字濾波器本質上是一個完成特定運算的數字計算過程，也可以理解為是一台計算機。描述離散系統輸出與輸入關系的卷積和差分方程只是給數字信號濾波器提供運算規則，使其按照這個規則完成對輸入數據的處理。時域離散系統的頻域特性:
其中
、
分別是數字濾波器的輸出序列和輸入序列的頻域特性（或稱為頻譜特性），
是數字濾波器的單位取樣響應的頻譜，又稱為數字濾經過濾波后的頻域響應。只要按照輸入信號頻譜的特點和處理信號的目的，適當選擇，使得濾波后滿足設計的要求，這就是數字濾波器的濾波原理。
數字濾波器根據其沖激響應函數的時域特性，可分為兩種，即無限長沖激響應(IIR)數字濾波器和有限長沖激響應(FIR)數字濾波器。IIR 數字濾波器的特征是，具有無限持續時間沖激響應，需要用遞歸模型來實現，其差分方程為：
系統函數為：
設計IIR濾波器的任務就是尋求一個物理上可實現的系統函數H(z)，使其頻率響應H(z)滿足所希望得到的頻域指標，即符合給定的通帶截止頻率、阻帶截止頻率、通帶衰減系數和阻帶衰減系數。
基本設計過程如下：
1.將給定的數字濾波器指標轉換成模擬濾波器的指標。
2.設計模擬濾波器。
3.將模擬濾波器轉換成數字濾波器系統函數。
（二）實驗步驟
1.在完成濾波器設計之后，采用filter()對輸入信號進行濾波處理。調用信號產生函數mstg產生由抑制載波調制信號相加構成的復合信號st，該函數還會自動繪圖顯示其時域波形和幅頻特性曲線，如下圖1所示。
由圖1中（a）和（b）可見，三路信號時域混疊無法在時域分離，但在頻域是可以分離的，可以通過濾波的方法進行分離，即通過設計IIR濾波器，分離這三個不同頻率的信號。
2.要求將三路信號進行分離，通過觀察st信號的幅頻特性曲線，分別確定可以分離st中三路抑制載波單頻調幅信號的三個濾波器（低通濾波器、帶通濾波器、高通濾波器）的通帶截止頻率和阻帶截止頻率。濾波器的通帶最大衰減為0.1dB，阻帶最小衰減為60dB。
3.調用ellipord()和ellip()分別設計這三個橢圓濾波器，並繪圖顯示其損耗函數曲線分別如圖2，圖3，圖4。
圖2 低通損耗函數曲線
4.調用filter()，用三個濾波器分別對信號產生函數mstg產生的信號進行濾波，分離出st中的不同載波頻率信號，並繪圖顯示三個信號的時域波形，分別如圖5，圖6，圖7。
三、實驗結果分析
（一）函數mstg
閱讀信號產生函數mstg，確定三路調幅信號的載波頻率和調制信號頻率。第1路調幅信號的載波頻率fc1=1000Hz，調制信號頻率fm1=100Hz。第2路調幅信號的載波頻率fc2=500Hz，調制信號頻率fm2=50Hz。第3路調幅信號的載波頻率fc2=250Hz，調制信號頻率fm2=25Hz。
- 采樣點數對頻譜圖的影響
1.調幅信號
當N=800,1600,1800,2000時，調幅信號產生的頻譜圖如圖8,9,10,11所示。
圖8 N=800時s（t）的波形及頻譜圖
分析：因為信號s（t）是周期序列，頻譜分析時要求觀察時間為整數倍周期。s（t）的每個頻率成分都是25Hz的整數倍。采樣頻率Fs=10kHz=25*400Hz，即在25Hz的正弦波的一個周期中采樣400個點。所以，當N為400的整數倍時一定為s（t）的整數倍周期。因此，采樣點數N=800,1600,2000時，對s（t）進行N點FFT可以得6根理想譜線，而當N=1800時，不是400的整數倍，則不能得到。
2.AM調幅信號
當N=800,1600,1800,2000時產生的頻譜圖如圖12,13,14,15所示。
分析：因為信號s（t）時周期序列，頻譜分析時要求觀察時間為整數倍周期。因此，采樣點數N=800,1600,2000時，對s（t）進行N點FFT可以得理想譜線，而當N=1800時，不是400的整數倍，則不能得到。當將該調幅信號修改為AM信號后，s（t）的頻譜中有較大的頻譜分量。如圖所示。
- IIR濾波器
- 濾波器參數選取
由（一）可知，三路調幅信號的載波頻率分別為250Hz，500Hz，1000Hz。帶寬為50Hz，100Hz，200Hz。所以分離混合信號st中三路抑制載波單頻調幅信號的三個濾波器（低通濾波器、帶通濾波器、高通濾波器）的指標參數選取如下：
對載波頻率為250Hz的調幅信號，用低通濾波器分離，其通帶截止頻率fp=280Hz，通帶最大衰減ap=0.1dB，阻帶截止頻率為fs=450Hz，阻帶最小衰減as=60dB。對載波頻率為500Hz的調幅信號，用低通濾波器分離，其通帶截止頻率fpl=440Hz、fph=560Hz，通帶最大衰減ap=0.1dB，阻帶截止頻率為fsl=275Hz、fsh=900Hz阻帶最小衰減as=60dB。對載波頻率為1000Hz的調幅信號，用低通濾波器分離，其通帶截止頻率fp=890Hz，通帶最大衰減ap=0.1dB，阻帶截止頻率為fs=600Hz，阻帶最小衰減as=60dB。
IIR.m
function main
st=mstg; %調用信號產生函數mstg產生由三路抑制載波調幅信號相加構成的復合信號s
Fs=10000; T=1/Fs; %采樣頻率
fp=280; fs=450;wp=2*fp/Fs;ws=2*fs/Fs; %低通
name='y_1(t)';LHL_filter(wp,ws,name,st,T,'low');
fp_L=440; fp_R=560;fs_L=275; fs_R=900;%帶通
wp=[2*fp_L/Fs,2*fp_R/Fs];ws=[2*fs_L/Fs,2*fs_R/Fs];
name='y_2(t)'; LHL_filter(wp,ws,name,st,T,'bandpass');
fp=890; fs=600;wp=2*fp/Fs;ws=2*fs/Fs; %高通
name='y_3(t)';LHL_filter(wp,ws,name,st,T,'high');
end
function LHL_filter(wp,ws,name,st,T,flag)
rp=0.1;rs=60; %DF指標（低通濾波器的通、阻帶邊界頻）
[N,wp]=ellipord(wp,ws,rp,rs); %調用ellipord計算橢圓DF階數N和通帶截止頻率wp?
[B,A]=ellip(N,rp,rs,wp,flag); %調用ellip計算橢圓帶通DF系統函數系數向量B和A
yt=filter(B,A,st); %濾波器軟件實現
figure; freqz(B,A);%B,A分別為系統函數分子，分母多項式系數向量?
figure;tplot(st,yt,T,name);
end
function tplot(st,xn,T,name)
%時域序列連續曲線繪圖函數?
%?xn:信號數據序列，%?T為采樣間隔
N=length(xn);n=0:N-1; t=n*T; Tp=N*T; f=n/Tp;
subplot(311);plot(t,xn);xlabel('t/s');ylabel(name);
axis([0,t(end),min(xn),1.2*max(xn)]);
subplot(312);stem(f,abs(fftshift(fft(st,N))));
title('原 st(t)的頻譜');xlabel('f/Hz');ylabel('幅度');
subplot(313);stem(f,abs(fftshift(fft(xn,N))));
title('濾波后xn(t)的頻譜');xlabel('f/Hz');ylabel('幅度');
end
function st=mstg
%產生信號序列向量st,並顯示st的時域波形和頻譜
%st=mstg 返回三路調幅信號相加形成的混合信號，長度N=1600
N=1600;Fs=10000;T=1/Fs;Tp=N*T; % N為信號st的長度。采樣頻率Fs=10kHz，Tp為采樣時間
t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;
fc1=Fs/10;fm1=fc1/10;%第1路調幅信號的載波頻率fc1=1000Hz,調制信號頻率fm1=100Hz
fc2=Fs/20;fm2=fc2/10;%第2路調幅信號的載波頻率fc2=500Hz,調制信號頻率fm2=50Hz
fc3=Fs/40;fm3=fc3/10;%第3路調幅信號的載波頻率fc3=250Hz,調制信號頻率fm3=25Hz
xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t); %產生第1路調幅信號
xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t); %產生第2路調幅信號
xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t); %產生第3路調幅信號
st=xt1+xt2+xt3; %三路調幅信號相加
fxt=fft(st,N); %計算信號st的頻譜
subplot(211);plot(t,st);xlabel('t/s');ylabel('s(t)');
axis([0,Tp/2,min(st),max(st)]);title('(a) s(t)的波形');
subplot(212);stem(f,abs(fxt)./max(abs(fxt)));
title('(頻譜');xlabel('f/Hz');ylabel('幅度');
endmstg.m
function st=mstg
%產生信號序列向量st,並顯示st的時域波形和頻譜
%st=mstg; 返回三路調幅信號相加形成的混合信號，長度N=1600
N=2000;Fs=10000;T=1/Fs;Tp=N*T; % N為信號st的長度。采樣頻率Fs=10kHz，Tp為采樣時間
t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;
fc1=Fs/10;fm1=fc1/10;%第1路調幅信號的載波頻率fc1=1000Hz,調制信號頻率fm1=100Hz
fc2=Fs/20;fm2=fc2/10;%第2路調幅信號的載波頻率fc2=500Hz,調制信號頻率fm2=50Hz
fc3=Fs/40;fm3=fc3/10;%第3路調幅信號的載波頻率fc3=250Hz,調制信號頻率fm3=25Hz
xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t); %產生第1路調幅信號
xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t); %產生第2路調幅信號
xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t); %產生第3路調幅信號
st=xt1+xt2+xt3; %三路調幅信號相加
fxt=fft(st,N); %計算信號st的頻譜
subplot(2,1,1);plot(t,st);xlabel('t/s');ylabel('s(t)');
title(['N=',num2str(N),'s(t)的波形']);axis([0,Tp/8,min(st),max(st)]);
subplot(2,1,2);stem(f,abs(fxt)/max(abs(fxt)));
title(['N=',num2str(N),'s(t)的頻譜']);xlabel('f/Hz');ylabel('幅度');
axis([0,Fs/5,0,1.2]);
endAM_mstg.m
function AM_mstg
N=2000; Fs=10000;T=1/Fs;Tp=N*T;t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;
fc1=Fs/10;fm1=fc1/10;%第1路調幅信號的載波頻率fc1=1000Hz,調制信號頻率fm1=100Hz
fc2=Fs/20;fm2=fc2/10;%第2路調幅信號的載波頻率fc2=500Hz,調制信號頻率fm2=50Hz
fc3=Fs/40;fm3=fc3/10;%第3路調幅信號的載波頻率fc3=250Hz,調制信號頻率fm3=25Hz
xt1=[1+cos(2*pi*fm1*t)].*cos(2*pi*fc1*t);
xt2=[1+cos(2*pi*fm2*t)].*cos(2*pi*fc2*t);
xt3=[1+cos(2*pi*fm3*t)].*cos(2*pi*fc3*t);
st=xt1+xt2+xt3;fxt=fft(st,N);
subplot(211);plot(t,st);xlabel('t/s');ylabel('s(t)');
axis([0,Tp/8,min(st),max(st)]);title(['N=',num2str(N),' AM s(t)波形圖'])
subplot(212);stem(f,abs(fxt)/max(abs(fxt)));title(['N= ',num2str(N),'頻譜圖'])
axis([0,Fs/5,0,1.2]);xlabel('f/Hz');ylabel('幅度');
end