FIR、IIR數(shù)字濾波器設(shè)計

1、實驗要求 1、給出音頻信號的時域和頻譜特性；2、設(shè)計一個IIR數(shù)字濾波器，給出濾波器的時域和頻譜特性，并利用濾波器對音頻信號進行濾波，給出濾波結(jié)果（濾波后的時域和頻譜特性）；3、設(shè)計一個FIR數(shù)字濾波器，給出濾波器的時域和頻譜特性，并利用濾波器對音頻信號進行濾波，給出濾波結(jié)果（濾波后的時域和頻譜特性）。實驗原理2.1 采樣頻率、位數(shù)及采樣定理 采樣頻率2，也稱為采樣速度或者采樣率，定義了每秒從連續(xù)信號中提取并組成離散信號的采樣個數(shù)，它用赫茲（Hz）來表示。采樣頻率的倒數(shù)是采樣周期或者叫作采樣時間，它是采樣之間的時間間隔。通俗的講采樣頻率是指計算機每秒鐘采集多少個聲音樣本，是描述聲音文件的音質(zhì)、

2、音調(diào)，衡量聲卡、聲音文件的質(zhì)量標(biāo)準(zhǔn)。采樣頻率越高，即采樣的間隔時間越短，則在單位時間內(nèi)計算機得到的聲音樣本數(shù)據(jù)就越多，對聲音波形的表示也越精確。 采樣位數(shù)可以理解為聲卡處理聲音的解析度。這個數(shù)值越大，解析度就越高，錄制和回放的聲音就越真實。我們首先要知道：電腦中的聲音文件是用數(shù)字0和1來表示的。所以在電腦上錄音的本質(zhì)就是把模擬聲音信號轉(zhuǎn)換成數(shù)字信號。反之，在播放時則是把數(shù)字信號還原成模擬聲音信號輸出。 采樣定理又稱奈奎斯特定理2，在進行模擬/數(shù)字信號的轉(zhuǎn)換過程中，當(dāng)采樣頻率fs不小于信號中最高頻率fm的2倍時，采樣之后的數(shù)字信號完整地保留了原始信號中的信息，一般實際應(yīng)用中保證采樣頻率為信號最高

3、頻率的510倍。 2.2時域信號的FFT分析 信號的頻譜分析就是計算機信號的傅里葉變換5。連續(xù)信號與系統(tǒng)的傅里葉分析顯然不便于用計算機進行計算，使其應(yīng)用受到限制。而FFT是一種時域和頻域均離散化的變換，適合數(shù)值運算，成為用計算機分析離散信號和系統(tǒng)的有力工具。對連續(xù)信號和系統(tǒng)，可以通過時域采樣，應(yīng)用DFT進行近似譜分析。 2.3 IIR數(shù)字濾波器設(shè)計原理 利用雙線性變換設(shè)計IIR濾波器（巴特沃斯數(shù)字低通濾波器的設(shè)計）6，首先要設(shè)計出滿足指標(biāo)要求的模擬濾波器的傳遞函數(shù)Ha(s),然后由Ha(s)通過雙線性變換可得所要設(shè)計的IIR濾波器的系統(tǒng)函數(shù)H(z)。如果給定的指標(biāo)為數(shù)字濾波器的指標(biāo)，則首先要轉(zhuǎn)

4、換成模擬濾波器的技術(shù)指標(biāo)，這里主要是邊界頻率Wp和Ws的轉(zhuǎn)換，對ap和as指標(biāo)不作變化。邊界頻率的轉(zhuǎn)換關(guān)系為 =2/T tan(w/2)。接著，按照模擬低通濾波器的技術(shù)指標(biāo)根據(jù)相應(yīng)設(shè)計公式求出濾波器的階數(shù)N和3dB截止頻率 c ；根據(jù)階數(shù)N查巴特沃斯歸一化低通濾波器參數(shù)表，得到歸一化傳輸函數(shù)Ha(p)；最后，將p=s/ c 代入Ha(p)去歸一，得到實際的模擬濾波器傳輸函數(shù)Ha(s)。之后，通過雙線性變換法轉(zhuǎn)換公式s=2/T(1-1/z)/(1+1/z)得到所要設(shè)計的IIR濾波器的系統(tǒng)函數(shù)H(z)。 2.4 FIR數(shù)字濾波器設(shè)計原理 基于窗函數(shù)的FIR 數(shù)字濾波器的設(shè)計方法通常也稱之為傅立葉級

5、數(shù)法，是用一定寬度窗函數(shù)截取無限脈沖響應(yīng)序列，獲得有限長的脈沖響應(yīng)序列，從而得到FIR 濾波器。它是在時域進行的，由理想濾波器的頻率響應(yīng))(wjdeH推導(dǎo)出其單位沖激響應(yīng)hd(n)，再設(shè)計一個FIR數(shù)字濾波器的單位沖激響應(yīng)h(n)去逼近hd(n)，表示此得到的離散濾波器的系統(tǒng)傳遞函數(shù)。該hd(n) 為無限長序列，因此Hd(z)是物理不可實現(xiàn)的。為了使系統(tǒng)變?yōu)槲锢砜蓪崿F(xiàn)的，且使實際的FIR 濾波器頻率響應(yīng)盡可能逼近理想濾波器的頻率響應(yīng)，采用窗函數(shù)將無限脈沖響應(yīng)hd(n)截取一段h(n)來近似表示hd(n)，可得：h (n) = hd(n)w(n) ,從而有：式中N 表示窗口長度，這樣H(z)就是

6、物理可實現(xiàn)的系統(tǒng)。并且從線性相位FIR濾波器的充要條件可知，為了獲得線性相位FIR 數(shù)字濾波器的沖激響應(yīng)h(n) ，那么序列h(n) 應(yīng)有= (N 1) / 2的延遲。由于窗函數(shù)的選擇對結(jié)果起著重要的作用，針對不同的信號和不同的處理目的來確定窗函數(shù)的選擇才能收到良好的效果。 實驗內(nèi)容1. 音頻信號的時域和頻譜特性x1=wavread(C:Documents and Settingsuser桌面新建文件夾chenwei.wav');figure(1);plot(x1);title(原始語音采樣后時域信號);xlabel(時間軸);ylabel(幅值);figure(2

7、);freqz(x1);title(原始語音信號采樣后頻譜);2. IIR數(shù)字濾波器fp = 4000; %通帶截止頻率fs = 4500; %阻帶截止頻率Fs=22050; %采樣點數(shù)wp=2*pi*fp/Fs; %通帶頻率歸一化ws=2*pi*fs/Fs; %阻帶頻率歸一化x1=wavread('C:Documents and Settingsuser桌面新建文件夾chenwei.wav');Rp=0.1; %通帶最大衰減Rs=50; %阻帶最小衰減Ts=1/Fs; %采樣頻率wp1=2/Ts*tan(wp/2); %將模擬指標(biāo)轉(zhuǎn)換成數(shù)字指標(biāo)ws1=2/

8、Ts*tan(ws/2); %將模擬指標(biāo)轉(zhuǎn)換成數(shù)字指標(biāo)N,Wn=buttord(wp1,ws1,Rp,Rs,'s'); %選擇濾波器的最小階數(shù)N以及濾波器固有頻率Wn %其中wp1和wp2以弧度為單位Z,P,K=buttap(N); %產(chǎn)生一個低通模擬butterworth濾波器原型Bap,Aap=zp2tf(Z,P,K); %零極點形式的傳遞函數(shù)的轉(zhuǎn)換b,a=lp2lp(Bap,Aap,Wn); %將表示為傳遞函數(shù)形式的模擬低通濾波器原型轉(zhuǎn) %成截止頻率為Wn的低通濾波器bd,ad=bilinear(b,a,Fs); %用雙線性變換法實現(xiàn)模擬濾波器到數(shù)字濾波器的轉(zhuǎn)換h,w=f

9、reqz(bd,ad); %數(shù)字濾波器的頻率響應(yīng)%畫濾波器特征曲線第9 / 10頁figure(1) subplot(1,1,1);plot(Fs*w/(2*pi),20*log10(abs(h); axis(0 8000 -150 20);grid;title('濾波器的性能分析');%畫濾波前后的時域語音信號figure(2)x,fs,bits=wavread('d:bobo.wav'); n=length(x);f=fs*(0:(n/2-1)/n;X=fft(x);z=filter(bd,ad,x);subplot(2,1,1);plot(x);title

10、('原始信號的波形');subplot(2,1,2);plot(z);title('濾波后信號的波形');%畫濾波前后的語音信號頻譜圖figure(3) sound(z,fs,bits);subplot(2,1,1);plot(f,abs(X(1:n/2);title('原始信號的頻譜');xlabel('Hz');Z=fft(z);subplot(2,1,2);plot(f,abs(Z(1:n/2);title('濾波后的信號頻譜')xlabel('Hz');分析 ： 可以看出，經(jīng)過濾波以后，高頻

11、部分已經(jīng)濾除。3. FIR數(shù)字濾波器clc;clear;close all;    fs = 8000;%采樣頻率為8000hz    x1=wavread('C:Documents and Settingsuser桌面新建文件夾chenwei.wav');%加入文件    t=0:1/8000:(size(x1)-1)/8000;    Au=0.05;%給語音信號加噪

12、聲    d=Au*os(2*pi*5000*t)'    x2=x1+d;    %濾波器的指標(biāo)    wp=0.25*pi;%通道邊界頻率    ws=0.3*pi;%阻帶截止頻率    Rp=1;%通帶最大衰減    Rs=15;%阻帶最小衰減    %利用窗函數(shù)設(shè)

13、計FIR濾波器    B = ws -  wp;    N = ceil(8*pi/B);    wc = (wp + ws) / 2 / pi;    b,a = fir1(N,wc,hamming(N+1);     

14、;   f2 = filter(b,a,x2);    figure(1);    freqz(b,a,512);    title('單位脈沖響應(yīng)的波形');    figure(2);    subplot(3,1,1);      plot(t,x1); %

15、畫出加躁前的時域圖    title('加躁前的時域波形');    subplot(3,1,2);    plot(t,x2); %畫出加躁后的時域圖    title('加躁后的時域波形');    subplot(3,1,3);    plot(t,f2); %畫出濾波后的時域圖  

16、;  title('濾波后的時域波形');    F0=fft(f2,1024);     f=fs*(0:511)/1024;    figure(3)      y2=fft(x2,1024);     subplot(2,1,1);      plot(f,abs(y2(1:512);      title('FIR低通濾波器濾波前的頻譜');     xlabel('頻率/Hz');     ylabel('幅值');     subp