2016年北航數字信號處理上機實驗二實驗報告.doc

離散時間信號處理實驗報告實驗二FIR數字濾波器的設計與實現班 級 學 號 姓 名 日 期 實驗二FIR數字濾波器的設計與實現實驗二FIR數字濾波器的設計與實現濾波器設計與實現是數字信號處理中最基本的內容之一。我們知道濾波器分IIR濾波器（無限沖激響應濾波器）和FIR濾波器（有限沖激響應濾波器）兩種。FIR濾波器的設計方法以直接逼近所需離散設計系統(tǒng)的頻率響應為基礎。FIR濾波器可以很容易地獲得線性相移特性，不存在不穩(wěn)定的問題，是實際系統(tǒng)中廣為采用的一種數字濾波器。FIR濾波器的設計，通常有窗函數設計法、頻率抽樣設計法和最佳逼近設計法。窗函數設計法比較簡單，它的頻率特性是理想濾波器頻譜與窗的頻譜的卷積，因而，其頻率特性取決于窗的類型和長度。頻率抽樣設計法比較直觀，但由于頻域的采樣會造成時域的混疊，從而濾波器疊性能不可能很高，為提高濾波器的性能，可以在過渡帶加上01之間的過渡點。本實驗主要采用窗函數設計法。FIR濾波器在Matlab中的實現包括時域卷積和filter濾波等方法，通過本實驗可以對兩種方法的實現過程和實現效果進行分析比較。相關Matlab函數說明1、有限沖激響應數字濾波器設計函數fir1和fir2Matlab函數fir1和fir2可以用來設計加窗的有限沖激響應數字濾波器。兩個函數均產生一個線性相位設計。函數fir1可用于設計常規(guī)的低通、高通、帶通和帶阻線性相位有限沖激響應濾波器。對于抽樣頻率為2Hz的情況，命令b = fir1(N,Wn)在向量b中返回以的升冪排列的N階低通或帶通濾波器的沖激響應系數。對于低通設計，歸一化截止頻率由標量Wn給定，它是在0和1之間的一個數。對于帶通設計，Wn是包含指定通帶邊界的一個雙元素向量Wn1,Wn2，其中0Wn1Wn21。命令b = fir1(N,Wn,high)其中N為一個偶數，用于設計高通濾波器。命令b = fir1(N,Wn,stop)其中Wn是一個雙元素向量，用于設計帶阻有限沖激響應濾波器。若沒有指定，用海明窗作為默認值。命令b = fir(N,Wn,taper)使用在向量taper中給定的長度為N+1的窗系數。然而，窗系數必須事先用一個適當的MATLAB函數來產生，如blackman，hamming，hanning，chebwin或kaiser。用到的命令為taper = blackman(N) taper = hamming(N) taper = hanning(N) taper = chebwin(N) taper = kaiser(N,beta)函數fir2可用于設計帶有任意形狀幅度響應的線性相位有限沖激響應濾波器。在其基本形狀中，命令是b = fir2(N,fpts,mval)它在長度為N+1的向量b中返回以的升冪排列的沖激響應系數。fpts是給定的頻率點向量，在0到1范圍內以遞增的順序排列，其中第一個頻率點是0，而最后的頻率點是1。和前面一樣，抽樣頻率被假定為2Hz。mval是給定頻率點的給定幅度值的向量，因此也必須和fpts的長度一樣，默認使用海明窗。為使用其他窗口，所用到的命令是b = fir2(N,fpts,mval,taper)其中向量taper包含給定的窗系數。2、FIR濾波器的實現(1) FIR濾波器的時域卷積實現濾波器是一種能夠讓某些頻率分量通過而拒絕其他頻率成分或對某些頻率成分進行修正的系統(tǒng)。對于一個線性時不變系統(tǒng)，其時域輸入和輸出關系可用卷積的方法描述，。對時域而言，均為數字信號，可在數字系統(tǒng)中實現。濾波器按照單位脈沖響應在某一有限區(qū)間外是否為零分為有限沖擊響應（FIR）濾波器和無限沖擊響應（IIR）濾波器。若單位脈沖響應在某一有限區(qū)間以外都是零，稱為有限沖擊響應FIR系統(tǒng)；若單位脈沖響應在某一有限區(qū)間外不是零，則稱為無限沖擊響應IIR系統(tǒng)。在線性時不變系統(tǒng)的卷積描述中，對于FIR濾波器，為有限項，范圍從0到M，故。即通過卷積運算，可實現對輸入信號的濾波處理。(2) FIR濾波器的filter濾波實現MATLAB中有濾波函數filter，用法為y=filter(b,a,x)。其中b和a為差分方程系數，x為輸入序列，y為濾波輸出序列。對于FIR濾波器，a設置為1即可。MATLAB實現主要包括兩個部分，即濾波器設計和濾波。代碼如下：%fir濾波器系數設計fl_kaiser = 750 1150;fl_mag = 1 0;fl_dev = 0.01 0.05;fl_n_kaiser,fl_wn,fl_beta,fl_ftype = kaiserord(fl_kaiser,fl_mag,fl_dev,fs);h = fir1(fl_n_kaiser,fl_wn,fl_ftype,kaiser(fl_n_kaiser+1,fl_beta);%濾波函數y = filter(h,1,x);3、wav格式音頻文件讀取函數 y = wavread(file)y,fs,nBits = wavread(file)輸入參數：file為wav格式的音頻文件。輸出參數：y為存放音頻數據的矩陣；fs為該音頻文件的采樣頻率；nBits表示語音信號的每個采樣點所占用的位數。4、wav格式音頻文件寫入函數wavwrite(y,fs,nBits,wavfile)wavwrite(y,fs,wavfile)wavwrite(y,wavfile)輸入參數：y是存放要寫入文件的音頻數據，fs是寫入文件的音頻信號采樣頻率，缺省的情況下默認fs為8000Hz。nBits是表示語音信號的每個取樣點所占用的位數，它必須是8,16,24或32。nBits缺省的情況下默認為16。wavfile是生成的WAV文件的文件名。實驗題目1、已知兩個信號序列：用conv函數求兩個序列的卷積和，并繪制三個序列的波形。（1）MATLAB 源代碼：clc;n1=0:20;n2=0:10;y1=0:30;x1n=0.8.n1;x2n=ones(1,11);y=conv(x1n,x2n);subplot(3,1,1)stem(n1,x1n);subplot(3,1,2)stem(n2,x2n);subplot(3,1,3)stem(y1,y);（2）實驗結果：第一張圖為f1，第二張圖為f2，第三張圖為f1*f22、利用wavread函數（或者audioread函數）將加噪聲的音頻文件“第一場雪加噪聲.wav”讀入Matlab的工作空間，利用load函數從h40.mat文件中讀取40階離散時間帶通濾波器的脈沖響應，用conv函數求兩個序列的卷積和，并利用wavwrite函數（或者audiowrite函數）將結果寫入音頻文件“第一場雪濾波40.wav”中。繪制兩個音頻信號的波形，并通過音頻播放器比較兩個音頻信號的區(qū)別。（1）MATLAB 源代碼：a,Fs,bits = wavread(第一場雪加噪聲);load(h40.mat);a=a;a1=a(1,1:end);a2=a(2,1:end);h1n=h40;Y1=conv(a1,h1n);subplot(6,1,1)plot(a1);xlabel(n1);ylabel(a1);subplot(6,1,2)plot(h1n);xlabel(n2);ylabel(h1n);subplot(6,1,3)plot(Y1);xlabel(y1);ylabel(a1*h1n);Y2=conv(a2,h1n);subplot(6,1,4)plot(a2);xlabel(n1);ylabel(a2);subplot(6,1,5)plot(h1n);xlabel(n2);ylabel(h1n);subplot(6,1,6)plot(Y2);xlabel(y2);ylabel(a2*h1n);ak=zeros(2,3109334);ak(1,1:end)=Y1;ak(2,1:end)=Y2;ak=ak;wavwrite(ak,Fs,bits,第一場雪濾波40.wav);（2）實驗結果：說明：n1為左聲道信號，n2為右聲道信號3、利用wavread函數（或者audioread函數）將加噪聲的音頻文件“第一場雪加噪聲.wav”讀入Matlab的工作空間，利用load函數從h392.mat文件中讀取392階離散時間帶通濾波器的脈沖響應，用conv函數求兩個序列的卷積和，并利用wavwrite函數（或者audiowrite函數）將結果寫入音頻文件“第一場雪濾波392.wav”中。繪制兩個音頻信號的波形，并通過音頻播放器比較兩個音頻信號的區(qū)別。（1）MATLAB 源代碼：a,Fs,bits = wavread(第一場雪加噪聲);load(h392.mat);a=a;a1=a(1,1:end);a2=a(2,1:end);n1=1:3109295;n2=1:393;y1=1:(3109295+393-1);h1n=h392;Y3=conv(a1,h1n);subplot(6,1,1)plot(n1,a1);xlabel(n1);ylabel(a1);subplot(6,1,2)plot(n2,h1n);xlabel(n2);ylabel(h1n);subplot(6,1,3)plot(y1,Y3);xlabel(y1);ylabel(a1*h1n);Y4=conv(a2,h1n);subplot(6,1,4)plot(n1,a2);xlabel(n1);ylabel(a2);subplot(6,1,5)plot(n2,h1n);xlabel(n2);ylabel(h1n);subplot(6,1,6)plot(y1,Y4);xlabel(y2);ylabel(a2*h1n);ak=zeros(2,3109687);ak(1,1:end)=Y3;ak(2,1:end)=Y4;ak=ak;wavwrite