一、實驗目的

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濾波器

1. 濾波器引數選取

由（一）可知，三路調幅訊號的載波頻率分別為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('幅度');
end

mstg.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]);
end

AM_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

資源點選下載