音樂信號濾波去噪-基于Gaussian的頻率采樣型FIR濾波器設(shè)計

1、 王杰 音樂信號濾波去噪-使用GAUSS窗設(shè)計的頻率采樣型FIR濾波器 第24頁 共23頁音樂信號濾波去噪基于Gausswin設(shè)計的FIR濾波器學(xué)生姓名：王杰 指導(dǎo)老師：高明摘 要 本次課程設(shè)計主要內(nèi)容是設(shè)計利用Gausssian窗口函數(shù)設(shè)計法設(shè)計一個FIR濾波器，將一段加入噪聲干擾的音樂信號進(jìn)行濾波去噪處理并根據(jù)濾波前后的波形和頻譜分析濾波性能。本課程設(shè)計仿真平臺為MATLAB7.0，開發(fā)工具是M語言編程。首先在因特網(wǎng)上下載工具一段音樂信號，并加入預(yù)設(shè)好頻率的單頻噪聲，然后對信號進(jìn)行頻譜分析得到加噪前后的信號對比圖，然后設(shè)計濾波器進(jìn)行濾波去噪處理，最后比較濾波前后的波形和頻譜并進(jìn)行分析。由分

2、析結(jié)果得到，濾波器后的語音信號與原始信號基本一致，即設(shè)計的FIR濾波器能夠去除信號中所加單頻噪聲，達(dá)到了設(shè)計目的。關(guān)鍵詞 MATLAB；濾波去噪；FIR濾波器；Gauss窗1引 言 本次課程設(shè)計主要是將一段音樂信號加入噪聲，然后用某種函數(shù)法設(shè)計出的FIR濾波器對加入噪聲后的音樂信號進(jìn)行濾波去噪處理，處理時采用的是利用窗口設(shè)計法選擇GAUSSWIN 設(shè)計的FIR 濾波器， 通過課程設(shè)計了解FIR 濾波器設(shè)計的原理和步驟，掌握用Matlab語言設(shè)計濾波器的方法，觀察音樂信號濾波前后的時域波形的比較，加深對濾波器作用的理解。1.1 課程設(shè)計的目的數(shù)字信號處理（Digital Signal Procc

3、essing，簡稱DSP）是一門涉及許多學(xué)科而廣泛應(yīng)用于許多領(lǐng)域的新型學(xué)科。20世紀(jì)60年代后，隨著計算機(jī)和信息技術(shù)的飛速發(fā)展，數(shù)字信號處理技術(shù)應(yīng)運而生并并得到迅速的發(fā)展。在過去的二十多年時間里，數(shù)字信號處理已經(jīng)在通信領(lǐng)域得到極為廣泛的應(yīng)用。數(shù)字信號處理是利用計算機(jī)或?qū)Ｓ锰幚碓O(shè)備，以數(shù)字形式對信號進(jìn)行采集、變換、濾波、估值、增強、壓縮、識別等處理，以得到符合人們的需要的信號形式。在本次課程設(shè)計中，最主要的設(shè)計是設(shè)計FIR 濾波器， FIR 濾波器的設(shè)計方法主要分為兩類，第一類是基于逼近理想濾波器特性的方法，包括窗函數(shù)法、頻域采樣法和等波紋最佳逼近法；第二類是最優(yōu)設(shè)計法。本次的課程設(shè)計主要采用的

4、是第一類設(shè)計方法，是利用GAUSSWIN法設(shè)計FIR濾波器對一段音樂信號進(jìn)行濾波去噪，通過這一過程，對濾波前后波形進(jìn)行對比分析得到結(jié)論。此課程設(shè)計比較簡單，主要是將書本中的知識應(yīng)用到現(xiàn)實中，并且根據(jù)自己對設(shè)計題目的理解，運用軟件編寫出程序?qū)崿F(xiàn)這一設(shè)計，也是我們對數(shù)字信號處理的原理進(jìn)行驗證的一個過程。對此，也可以加深我們對所學(xué)知識的理解，培養(yǎng)我們的動手能力。1.2 課程設(shè)計的要求（1）濾波器指標(biāo)必須符合工程實際。（2）設(shè)計完后應(yīng)檢查其頻率響應(yīng)曲線是否滿足指標(biāo)。（3）處理結(jié)果和分析結(jié)論應(yīng)該一致，而且應(yīng)符合理論。（4）獨立完成課程設(shè)計并按要求編寫課程設(shè)計報告書。1.3 設(shè)計平臺MATLABMATLA

5、B是由美國Math Works公司20世紀(jì)80年代中期推出的數(shù)學(xué)軟件。MATLAB是“Matric Laboratory”的縮寫，意及“矩陣實驗室”，優(yōu)秀的數(shù)值計算能力和卓越的數(shù)據(jù)可視化能力使其很快在數(shù)學(xué)軟件中脫穎而出。Matlab已經(jīng)發(fā)展成為多學(xué)科、多種工作平臺的功能強大的大型軟件。在歐美的高校和研究機(jī)構(gòu)中，MATLAB是一種非常流行的計算機(jī)語言，許多重要的學(xué)術(shù)刊物上發(fā)表的論文均是用MATLAB來分析計算以及繪制出各種圖形。MATLAB是一完整的并可擴(kuò)展的計算機(jī)環(huán)境，是一種進(jìn)行科學(xué)和工程計算的交互式程序語言。它的基本數(shù)據(jù)單元是不需要指定維數(shù)的矩陣，它可直接用于表達(dá)數(shù)學(xué)的算式和技術(shù)概念，而普通

6、的高級語言只能對一個個具體的數(shù)據(jù)單元進(jìn)行操作。因此，解決同樣的數(shù)值計算問題，使用MATLAB 要比使用Basic、Fortran 和C 語言等提高效率許多倍。許多人贊譽它為萬能的數(shù)學(xué)“演算紙”。MATLAB 采用開放式的環(huán)境，你可以讀到它的算法，并能改變當(dāng)前的函數(shù)或增添你自己編寫的函數(shù)MATLAB 包含的內(nèi)容非常豐富，功能強大，可以概括為以下幾個方面：(1)可以在多種操作系統(tǒng)下運行，如DOS、Windows 95/98/2000/2000/NT、CompaqAlpha、LinuxSun Solaris 等。(2)有超過500 種的數(shù)學(xué)、統(tǒng)計、科學(xué)及工程方面的函數(shù)，使用簡單快捷，并且有很強的用戶

7、自定義函數(shù)的能力。(3)有強大的圖形繪制和可視化功能， 可以進(jìn)行視覺數(shù)據(jù)處理和分析，進(jìn)行圖形、圖像的顯示及編輯，能夠繪制二維、三維圖形，使用戶可以制作高質(zhì)量的圖形，從而寫出圖文并茂的文章。(4)有從外部文件及外部硬件設(shè)備讀入數(shù)據(jù)的能力。(5)有豐富的工具箱toolbox。各個領(lǐng)域的專家學(xué)者將眾多學(xué)科領(lǐng)域中常用的算法編寫為一個個子程序，即m 文件，這些m 文件包含在一個個工具箱中。其工具箱可以分為兩大類，即功能性工具箱和科學(xué)性工具箱。功能性工具箱主要用來擴(kuò)充MATLAB的符號計算、圖形可視化、建模仿真、文字處理等功能以及與硬件實時交互的功能。學(xué)科性工具箱是按學(xué)科領(lǐng)域來分類的，如信號處理、控制、通

8、信、神經(jīng)網(wǎng)絡(luò)圖像處理、系統(tǒng)辨識、魯棒控制、模糊邏輯、小波等工具箱。MATLAB 中的信號處理工具箱內(nèi)容豐富，使用簡便。在數(shù)字信號處理中常用的算法，如FFT，卷積，相關(guān)，濾波器設(shè)計，參數(shù)模型等，幾乎都只用一條語句即可以調(diào)用。數(shù)字信號處理所常用的函數(shù)有波形的產(chǎn)生、濾波器的分析和設(shè)計、傅里葉變換、Z 變換等。2 設(shè)計原理2.1 FIR濾波器濾波器根據(jù)其沖激響應(yīng)函數(shù)的時域特性，可分為2 種，即無限長沖激響應(yīng)（IIR）濾波器和有限長沖激響應(yīng)（FIR）濾波器。FIR 和IIR 的濾波原理都是進(jìn)行卷積，就是對輸入信號進(jìn)行某種計算。FIR 用處就在于對數(shù)字信號進(jìn)行必要的處理，得到所需的輸出信號。FIR 系統(tǒng)有

9、自己突出的優(yōu)點：系統(tǒng)總是穩(wěn)定的；易實現(xiàn)線性相位；允許設(shè)計多通帶（或多阻帶）濾波器，后兩項都IIR系統(tǒng)不易實現(xiàn)的。FIR 數(shù)字濾波器的設(shè)計方法有多種，如窗函數(shù)設(shè)計法、頻率采樣法和Chebyshev逼近法等。隨著Matlab 軟件尤其是Matlab的信號處理工作箱的不斷完善，不僅數(shù)字濾波器的計算機(jī)輔助設(shè)計有了可能，而且還可以使設(shè)計達(dá)到最優(yōu)化。FIR 數(shù)字濾波器設(shè)計的基本步驟如下：（1）確定技術(shù)指標(biāo)在設(shè)計一個濾波器之前，必須首先根據(jù)工程實際的需要確定濾波器的技術(shù)指標(biāo)。在很多實際應(yīng)用中，數(shù)字濾波器常被用來實現(xiàn)選頻操作。因此，指標(biāo)的形式一般在頻域中給出幅度和相位響應(yīng)。幅度指標(biāo)主要以2種方式給出。第一種是

10、絕對指標(biāo)。他提供對幅度響應(yīng)函數(shù)的要求，一般應(yīng)用于FIR濾波器的設(shè)計。第二種指標(biāo)是相對指標(biāo)。他以分貝值的形式給出要求。本文中濾波器的設(shè)計就以線性相位FIR 濾波器的設(shè)計為例。（2）逼近確定了技術(shù)指標(biāo)后，就可以建立一個目標(biāo)的數(shù)字濾波器模型（通常采用理想的數(shù)字濾波器模型）。之后，利用數(shù)字濾波器的設(shè)計方法（窗函數(shù)法、頻率采樣法等） ，設(shè)計出一個實際濾波器模型來逼近給定的目標(biāo)。（3）性能分析和計算機(jī)仿真上兩步的結(jié)果是得到以差分或系統(tǒng)函數(shù)或沖激響應(yīng)描述的濾波器。根據(jù)這個描述就可以分析其頻率特性和相位特性，以驗證設(shè)計結(jié)果是否滿足指標(biāo)要求；或者利用計算機(jī)仿真實現(xiàn)設(shè)計的濾波器，再分析濾波結(jié)果來判斷。2.2 窗口

11、設(shè)計法窗口法設(shè)計的基本想法是要選取某一種合適的理想頻率選擇性濾波器（這種濾波器總是有一個非因果，無限長的脈沖響應(yīng)） ，然后將它的脈沖響應(yīng)截斷（或加窗）以得到一個線性相位和因果的FIR 濾波器。因此，這種方法的重點在于選擇某種恰當(dāng)?shù)拇昂瘮?shù)和一個合適的理想濾波器?，F(xiàn)用 代表一理想頻率選擇性濾波器，它在整個通帶內(nèi)有單位幅度增益和線性相位特性，而阻帶內(nèi)有零響應(yīng)。用窗口設(shè)計法基本步驟如下1：(1) 給定要求的理想頻率響應(yīng),一般給定分段常數(shù)的理想頻率特性。(2) 由于是在時域設(shè)計故必須求出 (2-1)(3) 由于 是無限時長的，故要用一個有限時長的“窗函數(shù)”序列 將加以截斷(相乘)，窗的點數(shù)是N點。截斷后

12、的序列為h(n)： (2-2)窗的點數(shù)N及窗的形狀是兩個極重要的參數(shù)。(4) 求出加窗后實際的頻率響應(yīng) (2-3)(5) 檢驗是否滿足的要求，不滿足，則需考慮改變窗形狀或改變 窗長的點數(shù)N，重復(fù)第(3)、(4)兩步，到滿足要求為止。常見的窗函數(shù)性能表如下表2-1所示。表2-1常見的窗函數(shù)性能表名稱濾波器過渡帶寬最小阻帶衰減名稱濾波器過渡帶寬最小阻帶衰減矩形1.8/M21dBPARZENWIN6.6/M56db巴特利特6.1/M25dBFLATTOPWIN19.6/M108db漢寧6.2/M44dBGAUSSWIN5.8/M30db漢明6.6/M51dBBARTHANNWIN3.6/M21db布

13、萊克曼11/M74dBBLACKMANHARRIS16.1/M109dbBOHMANWIN10/M48dbCHEBWIN15.2/M113dbNUTTALLWIN15.4/M108dbTUKEYWIN2.4/M22db2.3 GAUSSIAN窗高斯窗是一種指數(shù)窗。主瓣較寬，故而頻率分辨力低；無負(fù)的旁瓣，第一旁瓣衰減達(dá)一55dB。常被用來截短一些非周期信號，如指數(shù)衰減信號等。對于隨時間按指數(shù)衰減的函數(shù)，可采用指數(shù)窗來提高信噪比。Gausswin的時域表達(dá)式可以表示為： (2-4)其中k與 的取值范圍必須滿足22.4濾波器結(jié)構(gòu)本次課程設(shè)計的濾波器采用的是頻率采樣型結(jié)構(gòu)。(1) 頻率采樣型結(jié)構(gòu)的導(dǎo)入

14、若FIR DF的沖激響應(yīng)為有限長(N點)序列h(n)，則有：圖2-1 關(guān)于h(n)的推導(dǎo)所以，對h(n)可以利用DFT得到H(k),再利用內(nèi)插公式： (2-5)來表示系統(tǒng)函數(shù)。(2) 頻率采樣型濾波器結(jié)構(gòu)由式2-5得到FIR濾波器的另外一種結(jié)構(gòu)：頻率采樣型結(jié)構(gòu)。它是由兩部分級聯(lián)而成。 (2-6)其中：級聯(lián)中的第一部分為梳狀濾波器 ，第二部分由N個諧振器組成的諧振貴。 (2-7)(3) 頻率采樣型結(jié)構(gòu)流圖如圖2-2所示圖2-2 頻率采樣型結(jié)構(gòu)流程圖(4) 頻率采樣型結(jié)構(gòu)特點它的系數(shù)H(k)直接就是在 處的頻率響應(yīng)。因此，控制濾波器的頻率響應(yīng)是很直接的。但所有的相乘系數(shù)及H(k)都是復(fù)數(shù)，應(yīng)將它們

15、先化成二階的實數(shù)，這樣乘起來比較復(fù)雜，增加了乘法次數(shù)及儲存量；且所有的諧振器的極點都是在單位圓上，由 決定?？紤]到系數(shù)量化的影響，當(dāng)系數(shù)量化時，極點會移動，有些極點就不能被梳狀濾波器的零點鎖抵消。(零點由延時單元絕對，不受量化的影響)系統(tǒng)就會變得不穩(wěn)定1。3設(shè)計步驟3.1設(shè)計流程圖濾波器設(shè)計流程圖如圖3-1所示。圖3-1 濾波器設(shè)計流程圖3.2下載并截取音樂信號從網(wǎng)上下載一段音樂，從中截取一段格式為.wav 的語音信號，時間為23秒，并加入噪聲。在MALAB 平臺上，觀察原始語音信號與加入噪聲后的時域和頻譜圖。源程序如下所示：x,fs,bits=wavread('G:數(shù)字信號處理課程設(shè)

16、計_rEfunk.wav');%讀取音樂信號的數(shù)據(jù)sound(x,fs,bits);N=length(x); %計算信號x的長度t=0:1/fs:(N-1)/fs; %計算時間范圍，樣本除以采樣頻率x=x'fn=2436; %單頻噪聲功率y=x+0.1*sin(fn*2*pi*t); %加單頻噪聲X=abs(fft(x);Y=abs(fft(y); %對原始信號和加噪信號進(jìn)行fft變換X=X(1:N/2);Y=Y(1:N/2); %截取前半部分deltaf=fs/N; %計算頻譜的譜線間隔f=0:deltaf:fs/2-deltaf; %計算頻譜的頻率范圍原始信號與含噪信號的時

17、域圖和頻域圖如圖3-2所示，由時頻域分析圖可以看出，音樂信號能量主要集中在1000Hz以內(nèi)。圖3-2 原始信號與含噪信號的時頻域分析圖3.3濾波器設(shè)計截取好原始信號，接下來的工作是設(shè)計一個Gaussian窗濾波器。利用公式1： RP=-min(db(1:1:ceil(wp1/dww)+1,ceil(wp2/dww)+1:1:end) (3-1)AS=-max(db(wsd/dww+1:wsu/dww+1) (3-2)分別得到通帶波紋及通帶最小阻帶衰減，并輸入相應(yīng)濾波器的頻率參數(shù)，接著用freqz_m求得濾波器的響應(yīng)頻率特性。源程序如下所示:fp1=1836;fs1=2406;fs2=2476;

18、fp2=3036; %帶阻濾波器設(shè)置指標(biāo)Rp=6;As=28;fcd=(fp1+fs1)/2; %計算上下邊帶中心頻率fcu=(fp2+fp2)/2;df=min(fs1-fp1),(fp2-fs2); wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; % 將Hz為單位的模擬頻率換算為rad為單位的數(shù)字頻率wsd=fs1/fs*2*pi;wsu=fs2/fs*2*pi;M=ceil(5.8*pi/dw)+1; %計算Gauss窗設(shè)計該濾波器時需要的階數(shù)n=0:M-1; % 定義時間范圍gausswf=gausswin(M); %利用gausswin函

19、數(shù)產(chǎn)生一個Gaussian窗函數(shù)hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); %計算理想阻帶濾波器的脈沖響應(yīng)h_bs=gausswf'.*hd_bs; %窗口法計算實際濾波脈沖響應(yīng)db,mag,pha,grd,w=freqz_m(h_bs,1); %計算濾波器的頻率特性wp1=fp1/fs*2*pi;wp2=fp2/fs*2*pi;dww=2*pi/1000;RP=-min(db(1:1:ceil(wp1/dww)+1,ceil(wp2/dww)+1:1:end) %計算通帶最大衰減AS=-max(db(wsd/dww+1:

20、wsu/dww+1) %計算出阻帶最小衰減濾波器幅度響應(yīng)圖、相位響應(yīng)圖、脈沖響應(yīng)圖如圖3-3所示。由圖可知，濾波器的實際通帶最大衰減高于設(shè)定的通帶最大衰減，而實際的阻帶最小衰減低于設(shè)定的阻帶最小衰減，性能達(dá)標(biāo)。圖3-3濾波器性能分析圖3.4信號濾波處理設(shè)計好濾波器后，我們要對語音信號進(jìn)行濾波，對比濾波前和濾波后的音樂信號。源程序如下所示：y_fil=fftfilt(h_bs,y); %用設(shè)計好的濾波器對y進(jìn)行濾波Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:length(Y_fil)/2);% %計算頻譜取前一半sound (y_fil,fs,bits); %音樂信號回

21、放wavwrite(y_fil,fs,'d:funk_filted.wav'); %產(chǎn)生去噪后的音樂文件原始音樂信號、濾波后信號的時域、頻域分析圖，縱坐標(biāo)為對數(shù)坐標(biāo)的去噪前后的幅度頻譜圖對比圖如圖3-4所。圖3-4 音樂信號時域頻域分析圖3.5結(jié)果分析：通過觀察上圖，濾波后的語音信號發(fā)生了衰減，說明濾波器起到了濾波作用，同時通過頻譜對比，可以看出濾波器濾掉了一部分頻率范圍內(nèi)的信號。分別聽原始語音和濾波后的語音信號，發(fā)現(xiàn)濾波后的語音信號噪聲減小了，同時原始信號強度稍有減弱，基本達(dá)到了濾波的效果。3.6濾波器結(jié)構(gòu)設(shè)計本次課程設(shè)計采用的是頻率采樣型的濾波器結(jié)構(gòu)，將按照之前步驟設(shè)計出的

22、直接型濾波器通過tf2fs函數(shù)轉(zhuǎn)換成頻率采樣型1。源程序如下所示：function C,B,A=tf2fs(h)%直接型到頻率抽樣型的轉(zhuǎn)換%-%C,B,A=tf2fs(h)%C=包含各并聯(lián)部分增益的行向量%B=包含按行排列的分子系數(shù)矩陣%A=包含按行排列的分母系數(shù)矩陣%h=FIR濾波器單位沖激響應(yīng)向量%-N=length(h);H=fft(h,N);r=input('r=');magH=abs(H);phaH=angle(H)'%check even or odd Nif (N=2*floor(N/2) L=N/2-1; A1=1,-r,0;1,r,0; C1=real

23、(H(1),real(H(L+2);else L=(N-1)/2; A1=1,-r,0; C1=real(H(1);endk=1:L'%初始數(shù)組B和數(shù)組AB=zeros(L,2);A=zeros(L,3);%計算分母系數(shù)A(1:L,1)=1;A(1:L,2)=-2*r*cos(2*pi*k/N);A(1:L,3)=r.2A=A;A1;%計算分子系數(shù)B(1:L,1)=cos(phaH(2:L+1);B(1:L,2)=-cos(phaH(2:L+1)-(2*pi*k/N);%計算增益系數(shù)C=2*magH(2:L+1),C1'根據(jù)計算得出的C(各并聯(lián)部分增益的行向量),B(分子系數(shù)向

24、量),A(分母系數(shù)向量)畫出頻率采樣型濾波器結(jié)構(gòu)圖，如圖3-4所示，C,B,A參數(shù)如圖3-5所示。圖3-4 頻率采樣型濾波器結(jié)構(gòu)圖圖3-5 濾波器的C，B，A參數(shù)表4出現(xiàn)的問題及解決方法1、 下載得到的音樂文件為.mp3文件，用錄音機(jī)軟件轉(zhuǎn)換格式時無法識別，后來用Adobe Audition 3.0軟件將下載得到的音樂文件轉(zhuǎn)換成.wav文件就能順利轉(zhuǎn)換文件格式了。2、 在畫原始信號，加噪后的信號及濾波信號的時頻域分析圖時，發(fā)現(xiàn)畫出的圖形坐標(biāo)系上下限不一致，導(dǎo)致不能產(chǎn)生直觀的對比效果。經(jīng)過參考matlab的help文檔之后得知可以采用axis函數(shù)手動確定坐標(biāo)系上下限，直觀地畫出了對比圖。3、 在

25、生成加噪之后的音樂文件時，按照老師的提示使用wavwrite()函數(shù)生成，但是在括號里填寫好目標(biāo)存儲目錄之后發(fā)現(xiàn)matlab報錯，不能成功的生成。后來經(jīng)常參考matlab的help文件，得知wavwrite()函數(shù)需要填寫源數(shù)據(jù)作為參數(shù)，于是將加噪之后的音頻源數(shù)據(jù)y作為參數(shù)，填寫到wavwrite()成功生成了加了噪聲之后的文件。4、 在調(diào)用自編函數(shù)計算理想阻帶濾波器的脈沖響應(yīng)及計算濾波器的頻率特性時，發(fā)現(xiàn)ideal_lp函數(shù)以及freqz_m函數(shù)報錯，仔細(xì)檢查才發(fā)現(xiàn)沒有將已經(jīng)編譯好的m文件放入matlab的m文件夾內(nèi)，導(dǎo)致無法調(diào)用自編函數(shù)。將編譯好的函數(shù)m文件復(fù)制進(jìn)目標(biāo)文件夾后，順利地進(jìn)行了

26、計算。5、 設(shè)計好濾波器后，對濾波器的性能進(jìn)行驗證時，不知道該如何檢驗。后來按照老師所給的提示，通過公式計算出濾波器的實際Rp以及As與預(yù)設(shè)的Rp和As進(jìn)行比較并且在性能曲線圖中畫出作為參考的預(yù)設(shè)的Rp,As與實際的性能曲線相比較判斷性能是否達(dá)標(biāo)，經(jīng)對比發(fā)現(xiàn)濾波器符合設(shè)計標(biāo)準(zhǔn)，順利的完成了濾波器性能的檢驗6、 進(jìn)行的最后的濾波器結(jié)構(gòu)圖設(shè)計時，發(fā)現(xiàn)由tf2fs函數(shù)計算得出的C,B,A矩陣長度過長，達(dá)到了幾百個，將會使得結(jié)構(gòu)圖的繪制非常困難。經(jīng)過仔細(xì)調(diào)整通帶和阻帶截止頻率將高斯窗的階數(shù)設(shè)定在了43階，將C,B,A矩陣減少到了22行，簡化了濾波器的結(jié)構(gòu)圖，使得繪制順利的完成了。5結(jié)束語這次課程設(shè)計，

27、是運用上學(xué)期的所學(xué)數(shù)字信號處理理論知識的一次實踐，加深了對所學(xué)知識的理解，所以這次課程設(shè)計對我來說是收獲頗多的。本次的課程設(shè)計，我選到的課題主要內(nèi)容是通過數(shù)字信號處理技術(shù)實際處理信號，這里就是指對音樂信號的處理。作為存儲于計算機(jī)中的音樂信號，其本身就是離散化了的向量，我們只需將這些離散的量提取出來，就可以對其進(jìn)行處理了。本次課程設(shè)計的主要設(shè)計平臺是MATLAB，通過MATLAB 里調(diào)用幾個命令函數(shù)，很輕易的在實際化音樂與數(shù)字信號的理論之間搭了一座橋。課題的特色在于它將音樂看作了一個向量，于是音樂數(shù)字化了，則可以完全利用數(shù)字信號處理的知識來解決。我們可以像給一般信號做頻譜分析一樣，來給音樂信號做

28、頻譜分析，也可比較容易地用數(shù)字濾波器來對音樂進(jìn)行濾波處理。改變參數(shù)，理論結(jié)合實際，分析各參數(shù)對圖形的影響，從而加深對各個參數(shù)的理解。在完成這次課程設(shè)計過程中學(xué)到了許多東西，進(jìn)一步理解了濾波器設(shè)計方法和各參數(shù)意義，通過分析信號時域和頻域的關(guān)系等，加深了對濾波性能的理解，而且學(xué)會了使用Matlab 一些基本函數(shù)，增加了進(jìn)一步學(xué)習(xí)Matlab 軟件的興趣。同時，通過本次課程設(shè)計，鍛煉了我的動手能力，和提高了我分析問題，解決問題的能力。同學(xué)們在一起思考問題，通過親自動手設(shè)計一個FIR 濾波器，發(fā)現(xiàn)了許多潛在的問題，而這些問題是在平常被我們所忽略的，甚至認(rèn)為是不成問題的。讓我再次感受到我們應(yīng)當(dāng)把所學(xué)知識

29、和實踐相結(jié)合，才能夠提高自己的知識水平，取得更大的發(fā)展。參考文獻(xiàn)1程佩青.數(shù)字信號處理教程M.北京：清華大學(xué)出版社，2013:307,476,278-283.2yangpan011.幾種常見的窗函數(shù)以及其Matlab實現(xiàn)EB/OL. 3林愛英, 谷小青, 鄭寶周. 用頻率采樣法設(shè)計FIR濾波器J. 現(xiàn)代電子技術(shù), 2010, 33(17):85-87.附錄1：音樂信號濾波去噪設(shè)計源程序清單% 程序名稱： funk_filter.m% 程序功能：采用基于gausswin的窗口設(shè)計法，設(shè)計FIR 濾波器對含噪音樂信號進(jìn)行濾波去噪處理。% 程序作者：王杰% 最后修改日期： 2017-2-28x,fs

30、,bits=wavread('G:數(shù)字信號處理課程設(shè)計_rEfunk.wav'); %讀取音樂信號的數(shù)據(jù)sound(x,fs,bits);N=length(x); %計算信號x的長度t=0:1/fs:(N-1)/fs; %計算時間范圍，樣本除以采樣頻率x=x'fn=2436; %單頻噪聲功率y=x+0.1*sin(fn*2*pi*t); %加單頻噪聲X=abs(fft(x);Y=abs(fft(y); %對原始信號和加噪信號進(jìn)行fft變換X=X(1:N/2);Y=Y(1:N/2); %截取前半部分deltaf=fs/N; %計算頻譜的譜線間隔f=0:deltaf:fs/

31、2-deltaf; %計算頻譜的頻率范圍figure(1);subplot(2,2,1),plot(t,x);grid on;axis(0,2.5,-2,2);title('原始音樂信號');xlabel('時間（單位：s）');ylabel('幅度');subplot(2,2,2),plot(f,X);grid on;title('音樂信號幅度譜圖');xlabel('頻率（單位：Hz）');ylabel('幅度譜');subplot(2,2,3),plot(t,y);grid on;axis(0

32、,2.5,-2,2);title('加入噪聲干擾后的音樂信號');xlabel('時間（單位：s）');ylabel('幅度');subplot(2,2,4),plot(f,Y);grid on;title('加入噪聲干擾后的音樂信號幅度譜圖');xlabel('頻率（單位：Hz）');ylabel('幅度譜');wavwrite(y,fs,'d:funk_converted.wav');fp1=1836;fs1=2406;fs2=2476;fp2=3036; %帶阻濾波器設(shè)置指標(biāo)R

33、p=6;As=28;fcd=(fp1+fs1)/2; %計算上下邊帶中心頻率fcu=(fp2+fp2)/2;df=min(fs1-fp1),(fp2-fs2); wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; %將Hz為單位的模擬頻率換算為rad為單位的數(shù)字頻率wsd=fs1/fs*2*pi;wsu=fs2/fs*2*pi;M=ceil(5.8*pi/dw)+1; %計算Gauss窗設(shè)計該濾波器時需要的階數(shù)n=0:M-1; % 定義時間范圍gausswf=gausswin(M); %利用gausswin函數(shù)產(chǎn)生一個Gaussian窗函數(shù)hd_bs=

34、ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); %計算理想阻帶濾波器的脈沖響應(yīng)h_bs=gausswf'.*hd_bs; %窗口法計算實際濾波脈沖響應(yīng)db,mag,pha,grd,w=freqz_m(h_bs,1); %計算濾波器的頻率特性wp1=fp1/fs*2*pi;wp2=fp2/fs*2*pi;dww=2*pi/1000;RP=-min(db(1:1:ceil(wp1/dww)+1,ceil(wp2/dww)+1:1:end) %計算通帶最大衰減AS=-max(db(wsd/dww+1:wsu/dww+1) %計算出阻帶最小衰減f

35、igure(2);subplot(2,2,1);plot(w/pi,db);xlabel('w/pi');ylabel('db');title('FIR濾波器的幅度響應(yīng)圖');axis(0.3 0.9 -80 5);line(0,1.5,-RP,-RP,'Color','r','LineWidth',2,'LineStyle','-'); line(0,1.5,-AS,-AS,'Color','r','LineWidth'

36、;,2,'LineStyle','-'); line(wp1/pi,wp1/pi,-100,10,'Color','r','LineWidth',2,'LineStyle','-'); line(wp2/pi,wp2/pi,-100,10,'Color','r','LineWidth',2,'LineStyle','-');line(wsd/pi,wsd/pi,-100,10,'Color'

37、;,'r','LineWidth',2,'LineStyle','-'); line(wsu/pi,wsu/pi,-100,10,'Color','r','LineWidth',2,'LineStyle','-'); grid on; subplot(2,2,2);plot(w/pi,mag);xlabel('w/pi');ylabel('幅度mag');title('FIR濾波器的幅度響應(yīng)圖');axi

38、s(0 1 -0.2 1.2);grid on; subplot(2,2,3);plot(w/pi,pha);xlabel('w/pi');ylabel('相位pha');title('濾波器相位響應(yīng)圖');axis(0 1 -4 4);grid on;subplot(2,2,4);stem(n,h_bs);xlabel('n');ylabel('h(n)');title('濾波器脈沖響應(yīng)圖');grid on;y_fil=fftfilt(h_bs,y); %用設(shè)計好的濾波器對y進(jìn)行濾波 Y_fil

39、=abs(fft(y_fil);Y_fil=Y_fil(1:length(Y_fil)/2); %計算頻譜取前一半 sound (y_fil,fs,bits); %音樂信號回放wavwrite(y_fil,fs,'d:funk_filted.wav'); %產(chǎn)生去噪后的音樂文件figure(3); %作圖subplot(4,2,1);plot(t,x);grid on; title('原始音樂信號');xlabel('時間（t）');ylabel('幅度');axis(0,2.5,-2,2);subplot(4,2,2);plot(f,X);grid on; axis(0,3000,0,1500);title('原始音樂信號幅度譜');xlabel('頻率（f）');ylabel('幅度譜'); subplot(4,2,5);plot(t,y_fil);grid on; axis(0,2.5,-2,2);title('濾波后的音樂信號');xlabel('時間（t）');ylabel('幅度');subplot(4,2,6);plot