音樂濾波去噪——使用flattopwin窗設(shè)計(jì)FIR濾波器

1、 蔣敏 音樂信號(hào)去噪使用FLATTOPWIN窗設(shè)計(jì)的FIR濾波器 第28頁，共28頁音樂信號(hào)濾波去噪使用FLATTOPWIN窗設(shè)計(jì)的FIR濾波器學(xué)生姓名： 指導(dǎo)老師：黃紅兵摘 要 本次課程設(shè)計(jì)是使用Flattopwin窗設(shè)計(jì)FIR濾波器對(duì)音樂信號(hào)進(jìn)行濾波去噪。通過MATLAB軟件，運(yùn)用窗函數(shù)法來設(shè)計(jì)濾波器。從網(wǎng)上下載一段滿足要求的音樂，為它加入噪聲信號(hào)，觀察加噪前后的頻譜，采用窗函數(shù)設(shè)計(jì)法，給定相應(yīng)的技術(shù)指標(biāo)，設(shè)計(jì)一個(gè)滿足要求的濾波器，對(duì)音樂信號(hào)進(jìn)行濾波去噪處理。比較原始音樂信號(hào)與濾波后的時(shí)域波形圖，頻譜圖，回放濾波后的音樂信號(hào)，可聽見濾波后的音樂信號(hào)與原始音樂信號(hào)無大致差別，成功的實(shí)現(xiàn)了濾波

2、達(dá)到了設(shè)計(jì)要求。關(guān)鍵詞 MATLAB；濾波去噪；FIR濾波器；Flattopwin窗； 1 引 言本次課程設(shè)計(jì)是通過窗函數(shù)法設(shè)計(jì)一個(gè)Flattopin的FIR濾波器對(duì)有噪聲的音樂信號(hào)進(jìn)行濾波去噪。在課程設(shè)計(jì)中，了解FIR濾波器的性能與原理，也了解他的設(shè)計(jì)方法和步驟。掌握了用MATLAB語言設(shè)計(jì)濾波器，通過觀察音樂信號(hào)濾波前后的時(shí)域信號(hào)以及頻譜更加具體的了解了濾波器的作用。1.1 課程設(shè)計(jì)目的通過利用MATLAB 軟件來利用Flattopin設(shè)計(jì)FIR濾波器對(duì)音樂信號(hào)濾波去噪。使得我們更加熟悉MATLAB的語言環(huán)境，更加熟悉MATLAB語言的編程規(guī)則。并且在課程設(shè)計(jì)中通過觀察濾波器的幅度，相位圖

3、對(duì)Flattowin有了更加深刻地了解。也在窗函數(shù)的設(shè)計(jì)過程中，對(duì)濾波器的性能，功能以及設(shè)計(jì)方法有著更具體的了解和體驗(yàn)。通過本次課程設(shè)計(jì)，增強(qiáng)了我們獨(dú)立解決問題的能力，提高了自己的動(dòng)手能力。1.2 課程設(shè)計(jì)要求從網(wǎng)上下載一段.wav格式的音樂，繪制觀察時(shí)域波形及頻譜圖。對(duì)音樂信號(hào)加入噪聲干擾，根據(jù)Flattopwin的性能指標(biāo)合理設(shè)計(jì)FIR濾波器，再用濾波器對(duì)加入干擾的音樂信號(hào)進(jìn)行濾波去噪，比較濾波前后的頻譜圖并進(jìn)行分析。再回放語音信號(hào)對(duì)比原語音信號(hào)，查看濾波器是否對(duì)語音信號(hào)進(jìn)行了濾波去噪。1.3課程設(shè)計(jì)平臺(tái)本次課程設(shè)計(jì)通過MATLAB實(shí)現(xiàn)，MATLA是美國MathWorks公司出品的商業(yè)數(shù)學(xué)

4、軟件，用于算法開發(fā)、數(shù)據(jù)可視化、數(shù)據(jù)分析以及數(shù)值計(jì)算的高級(jí)技術(shù)計(jì)算語言和交互式環(huán)境，主要包括MATLAB和Simulink兩大部分。MATLAB是matrix&laboratory兩個(gè)詞的組合，意為矩陣工廠（矩陣實(shí)驗(yàn)室）。是由美國mathworks公司發(fā)布的主要面對(duì)科學(xué)計(jì)算、可視化以及交互式程序設(shè)計(jì)的高科技計(jì)算環(huán)境。它將數(shù)值分析、矩陣計(jì)算、科學(xué)數(shù)據(jù)可視化以及非線性動(dòng)態(tài)系統(tǒng)的建模和仿真等諸多強(qiáng)大功能集成在一個(gè)易于使用的視窗環(huán)境中，為科學(xué)研究、工程設(shè)計(jì)以及必須進(jìn)行有效數(shù)值計(jì)算的眾多科學(xué)領(lǐng)域提供了一種全面的解決方案，并在很大程度上擺脫了傳統(tǒng)非交互式程序設(shè)計(jì)語言（如C、Fortran）的編輯模

5、式，代表了當(dāng)今國際科學(xué)計(jì)算軟件的先進(jìn)水平1。MATLAB軟件包括五大通用功能：數(shù)值計(jì)算功能（Nemeric）；符號(hào)運(yùn)算功能（Symbolic）；數(shù)據(jù)可視化功能（Graphic）；數(shù)據(jù)圖形文字統(tǒng)一處理功能（Notebook）和建模仿真可視化功能（Simulink）。其中，符號(hào)運(yùn)算功能的實(shí)現(xiàn)是通過請(qǐng)求MAPLE內(nèi)核計(jì)算并將結(jié)果返回到MATLAB命令窗口。該軟件有三大特點(diǎn)：一是功能強(qiáng)大；二是界面友善、語言自然；三是開放性強(qiáng)。MATLAB在線性代數(shù)、矩陣分析、數(shù)值及優(yōu)化、數(shù)理統(tǒng)計(jì)和隨機(jī)信號(hào)分析、電路與系統(tǒng)、系統(tǒng)動(dòng)力學(xué)、信號(hào)和圖像處理、控制理論分析和系統(tǒng)設(shè)計(jì)、過程控制、建模和仿真、通信系統(tǒng)、以及財(cái)政金融

6、等眾多領(lǐng)域的理論研究和工程設(shè)計(jì)中得到了廣泛應(yīng)用。MATLAB在信號(hào)與系統(tǒng)中的應(yīng)用主要包括符號(hào)運(yùn)算和數(shù)值計(jì)算仿真分析。由于信號(hào)與系統(tǒng)課程的許多內(nèi)容都是基于公式演算，而MATLAB借助符號(hào)數(shù)學(xué)工具箱提供的符號(hào)運(yùn)算功能能基本滿足信號(hào)與系統(tǒng)課程的需求。例如，解微分方程、傅里葉正反變換、拉普拉斯正反變換、z正反變換等。MATLAB在信號(hào)與系統(tǒng)中的另一主要應(yīng)用是數(shù)值計(jì)算與仿真分析，主要包括函數(shù)波形繪函數(shù)運(yùn)算、沖激響應(yīng)與階躍響應(yīng)仿真分析、信號(hào)的時(shí)域分析、信號(hào)的頻譜分析、系統(tǒng)的S域分析、零極點(diǎn)圖繪制等內(nèi)容。數(shù)值計(jì)算仿真分析可以幫助學(xué)生更深入理解信號(hào)與系統(tǒng)的理論知識(shí)，并為將來使用MATLAB進(jìn)行信號(hào)處理領(lǐng)域的各

7、種分析和實(shí)際應(yīng)用打下基礎(chǔ)2。2 基本理論2.1 FIR濾波器FIR濾波器：有限長單位沖激響應(yīng)濾波器，是數(shù)字信號(hào)系統(tǒng)中最基本的元件，它可以在保證任意幅頻特性的同時(shí)具有嚴(yán)格的線性相頻特性，同時(shí)其單位抽樣響應(yīng)是有限長的，因而濾波器是穩(wěn)定的系統(tǒng)。因此，F(xiàn)IR濾波器在通信、圖像處理、模式識(shí)別等領(lǐng)域都有著廣泛的應(yīng)用。FIR濾波器的基本結(jié)構(gòu)：橫截型，級(jí)聯(lián)型，頻率抽樣型，快速卷積型，在硬件上一般通過集成電路，DSP芯片，可編程邏輯器件，F(xiàn)PGA/CPLD來實(shí)現(xiàn)3。FIR數(shù)字濾波器的特點(diǎn)：優(yōu)點(diǎn)：(1) 很容易獲得嚴(yán)格的線性相位，避免所處理的信號(hào)產(chǎn)生相位失真； (2) 極點(diǎn)全部在原點(diǎn)，無穩(wěn)定性問題； (3) 任何

8、一個(gè)非因果的有限長序列，總可以通過一定的延時(shí)，轉(zhuǎn)化為因果序列，所以因果性總是可以滿足； (4) 無反饋運(yùn)算，運(yùn)算誤差小。缺點(diǎn)：(1) 因?yàn)闊o極點(diǎn)，要獲得好的過渡帶特性，要以較高的結(jié)束為代價(jià)； (2) 無法利用模擬濾波器的設(shè)計(jì)結(jié)果，一般無解析設(shè)計(jì)公式，要借助輔助程序設(shè)計(jì)完成。一般來說濾波器的設(shè)計(jì)在于尋找一個(gè)頻率響應(yīng)去逼近，逼近的方法有三種窗函數(shù)設(shè)計(jì)法（時(shí)域逼近）；頻率采樣設(shè)計(jì)法(頻域逼近)；最優(yōu)化設(shè)計(jì)法（等波紋逼近）。本次課程設(shè)計(jì)采用窗函數(shù)設(shè)計(jì)法4。2.2 窗口設(shè)計(jì)法窗函數(shù)設(shè)計(jì)法是從單位脈沖響應(yīng)序列著手，使得逼近理想單位脈沖響應(yīng)序列 (2.1)窗函數(shù)設(shè)計(jì)FIR濾波器的步驟：(1) 根據(jù)技術(shù)要求確

9、定濾波器的頻率響應(yīng)確定其對(duì)應(yīng)的單位脈沖響應(yīng)。 (2.2)(2) 根據(jù)對(duì)過渡帶及阻帶衰減指標(biāo)的要求，選擇窗函數(shù)形式，并估計(jì)窗口長度N。(3) 計(jì)算濾波器的單位取樣響應(yīng)。 (2.2)(4) 計(jì)算濾波器的頻率響應(yīng)。 (2.3) 改變窗函數(shù)的形狀，可改善濾波器的特性，窗函數(shù)有許多種，但要滿足以下兩點(diǎn)： (1) 窗譜主瓣寬度要窄，已獲得較陡的過渡帶； (2) 相對(duì)于主瓣寬度，旁瓣要盡可能小，是能量盡量集中在主瓣中，這樣就可以減小肩峰和余振，以提高阻帶衰減和通帶平穩(wěn)性。 但實(shí)際上這兩點(diǎn)不能兼得，一般總是通過增加主瓣寬度來換取對(duì)旁瓣的抑制5。常見的窗函數(shù)性能表如表2.1所示：表2.1常見窗函數(shù)性能表名稱濾波

10、器過渡帶寬最小阻帶衰減名稱濾波器過渡帶寬最小阻帶衰減矩形1.8/M21dBPARZENWIN6.6/M56dB巴特利特6.1/M25dBFLATTOPWIN19.6/M108dB漢寧6.2/M44dBGAUSSWIN5.8/M60dB漢明6.6/M51dBBARTHANNWIN3.6/M40dB布萊克曼11/M74dBBLACKMANHARRIS16.1/M109dBBOHMANWIN5.8/M51.5dBCHEBWIN15.2/M113dBNUTTALLWIN15.4/M108dBTUKEYWIN2.4/M22dB2.3 FLATTOPWIN窗w=Flattopwin (L) 返回L-點(diǎn)Fl

11、attopwin窗口中列向量。Flattopwin窗的濾波器的過渡帶寬為19.6/M，最小阻帶衰減108db。定義式： (2.4)時(shí)間波形和幅度譜：圖2.1 Flattopwin窗的時(shí)間波形圖2.2 Flattopwin窗的幅度譜2.4濾波器的結(jié)構(gòu)(1) 橫截型差分方程： (2.5)圖2.3 FIR濾波器橫截型結(jié)構(gòu)(2)級(jí)聯(lián)型將差分方程分解為實(shí)系數(shù)二階因式乘積形式： (2.6)圖2.4 FIR濾波器級(jí)聯(lián)型結(jié)構(gòu)當(dāng)N為偶數(shù)時(shí)，其中有一個(gè)=0（N-1個(gè)零點(diǎn)）。級(jí)聯(lián)型每個(gè)基本節(jié)控制一對(duì)零點(diǎn)，便于控制濾波器的傳輸零點(diǎn)；系數(shù)比直接型多，所需乘法運(yùn)算多。(3)線性相位型FIR 的重要特點(diǎn)是可設(shè)計(jì)成具有嚴(yán)格線

12、性相位的濾波器，此時(shí)滿足偶對(duì)稱或奇對(duì)稱。 偶對(duì)稱時(shí)： (2.7)N為奇對(duì)稱時(shí)： (2.8)圖2.5 N為偶數(shù)時(shí)的線性相位結(jié)構(gòu)FIR濾波器圖2.6 奇為偶數(shù)時(shí)的線性相位結(jié)構(gòu)FIR濾波器(5) 頻率采樣型H(z)由兩部分組成：第一部分（FIR部分）： (2.9) 第二部分（IIR部分）： (2.10)圖2.7 FIR濾波器的頻率采樣型結(jié)構(gòu)這一結(jié)構(gòu)的最大的特點(diǎn)是它的系數(shù)H(k)直接就是濾波器在處的響應(yīng)，因此控制濾波器響應(yīng)很直接6。 3 設(shè)計(jì)步驟3.1 設(shè)計(jì)流程圖根據(jù)設(shè)計(jì)的要求，從網(wǎng)上下載一段音樂信號(hào)，在通過系統(tǒng)自帶的軟件對(duì)音樂信號(hào)進(jìn)行格式轉(zhuǎn)換以及剪裁，然后對(duì)音樂信號(hào)加入噪聲干擾，再利用Flattop

13、win窗設(shè)計(jì)合理的FIR濾波器。最后用濾波器對(duì)干擾后的音樂信號(hào)進(jìn)行濾波去噪。具體設(shè)計(jì)流程圖如下圖3.1所示：圖3.1設(shè)計(jì)流程圖 3.2 錄制音樂信號(hào)從電腦上下載一段格式為.wav的純音樂，并命名為yinyue.wav，通過錄音機(jī)將原始音樂信號(hào)的格式更改為PCM編碼，采樣率為8KHz，8位，單聲道，7KB/s并將其裁剪為10s的文件，調(diào)用wavread函數(shù)讀取信號(hào)的參數(shù)，繪制出音樂信號(hào)的時(shí)域波形。 圖3.2用錄音機(jī)將音樂進(jìn)行裁剪 圖3.3用錄音機(jī)將音樂格式進(jìn)行轉(zhuǎn)化圖3.4 原始音樂信號(hào)時(shí)域波形圖繪制音樂信號(hào)的時(shí)域波形后，對(duì)音樂信號(hào)進(jìn)行快速傅里葉變換，得到信號(hào)的頻譜特性，繪制出音樂信號(hào)的頻譜圖，向

14、音樂信號(hào)加入噪聲，并對(duì)加入噪聲的信號(hào)進(jìn)行快速傅里葉變換，繪制出加入噪聲后的時(shí)域波形，以及頻譜圖，將加入噪聲后的信號(hào)時(shí)域波形圖以及頻譜圖與未加入之前的波形圖和頻譜圖進(jìn)行對(duì)比。>>x=x' y=x+0.1*sin(fn*2*pi*t); %在原信號(hào)中加入噪聲>>X=abs(fft(x); Y=abs(fft(y); % 對(duì)原始信號(hào)和加噪信號(hào)進(jìn)行fft變換,取幅度譜>>X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分>>deltaf=fs/N; % 計(jì)算頻譜的譜線間隔>>f=0:deltaf:fs/2-deltaf;

15、 % 計(jì)算頻譜頻率范圍圖3.5 干擾前后信號(hào)時(shí)域圖與頻譜圖的比較3.3 濾波器設(shè)計(jì) 根據(jù)加入噪聲干擾后的頻譜圖，給出濾波器的性能指標(biāo)，在通過窗函數(shù)法設(shè)計(jì)得到自己的濾波器，得到數(shù)字濾波器的參數(shù)a，b。其中a為系統(tǒng)函數(shù)的分母系數(shù)，b為系統(tǒng)函數(shù)的分子系數(shù)。在通過調(diào)用函數(shù)即可求得濾波器的頻率響應(yīng)。在計(jì)算濾波器的實(shí)際AS以及RP查看濾波器是否滿足設(shè)計(jì)指標(biāo)。>>fpd=1825;fsd=2075;fsu=2325;fpu=2575;Rp=1;As=20; % 帶阻濾波器設(shè)計(jì)指標(biāo)>>RP=-min(db(1:wpd/dw+1,ceil(wpu/dw)+1:end) %計(jì)算濾波器的指標(biāo)

16、檢查是否達(dá)標(biāo)>>AS=-max(db(wsd/dw+1:wsu/dw+1)通過計(jì)算可以得出RP=0.7939，AS=20.3545滿足指標(biāo)要求，則濾波器符合要求。圖3.6設(shè)計(jì)的濾波器圖形3.4 信號(hào)濾波處理 濾波器設(shè)計(jì)完成后，在MATLAB平臺(tái)上用函數(shù)filter對(duì)加有噪聲的音樂信號(hào)進(jìn)行濾波，并求得濾波后的頻譜。繪制濾波前后音樂信號(hào)的時(shí)域波形圖和頻譜圖對(duì)比圖，以及原信號(hào)的時(shí)域波形圖及頻譜圖。 >>y_fil=filter(h_bs,1,y);%用設(shè)計(jì)好的濾波器對(duì)y進(jìn)行濾波圖3.7 濾波前后音樂信號(hào)的時(shí)域波形圖和頻譜圖對(duì)比圖3.5 結(jié)果分析 觀察濾波前后的頻譜圖，可明顯看

17、看到干擾噪聲的頻譜被濾除，觀察濾波后的時(shí)域波形圖與原信號(hào)的波形圖基本一致，在通過sound()函數(shù)對(duì)經(jīng)過Flattopwin窗設(shè)計(jì)的FIR濾波器之后的音樂信號(hào)進(jìn)行回訪可以聽出濾波之后的信號(hào)跟原始信號(hào)一樣清晰，完全濾除掉了噪聲的干擾。>> sound (y_fil,fs,bits); %播放濾波后的音樂信號(hào)>>sound(x,fs,bits); % 按指定的采樣率和每樣本編碼位數(shù)回放所得結(jié)果證明了用Flattopwin窗設(shè)計(jì)的FIR濾波器對(duì)音樂信號(hào)進(jìn)行濾波去噪是成功的。3.6 濾波器結(jié)構(gòu)設(shè)計(jì)在繪制濾波器結(jié)構(gòu)時(shí)由于我選擇的級(jí)聯(lián)型，因此要計(jì)算相關(guān)參數(shù)，可以通過相關(guān)函數(shù)來進(jìn)行計(jì)

18、算>>C,B,A=DIR2CAS(h_bs,1); %計(jì)算級(jí)聯(lián)型濾波器的參數(shù)然后在Viso中繪制出FIR濾波器的級(jí)聯(lián)型結(jié)構(gòu)即可。圖3.6 FIR濾波器的級(jí)聯(lián)型結(jié)構(gòu)4出現(xiàn)的問題及解決方法在本次課程設(shè)計(jì)中我遇到的問題如下：1、 繪制加入噪聲干擾的頻譜圖是不知如何調(diào)控縱坐標(biāo)，導(dǎo)致頻譜圖十分的不理想2、 在達(dá)到設(shè)計(jì)指標(biāo)且完美的濾除噪聲干擾信號(hào)時(shí)，濾波器的階數(shù)總是過大，而將階數(shù)降低時(shí)，濾波器總是不能比較完美的濾除噪聲3、 繪制濾波器結(jié)構(gòu)時(shí)選擇的是級(jí)聯(lián)型，但是卻找不到，計(jì)算級(jí)聯(lián)型參數(shù)的函數(shù)針對(duì)以上問題，相應(yīng)的解決方案如下：1、 通過向同學(xué)詢問了解到有函數(shù)可以用來調(diào)控橫縱坐標(biāo)2、 通過不斷設(shè)置濾

19、波器參數(shù)值最終設(shè)計(jì)出理想的濾波器，可在階數(shù)較低的情況下完美的濾出噪聲，也達(dá)到了設(shè)計(jì)指標(biāo)。3、 在查閱歷屆的PPT時(shí)，發(fā)現(xiàn)了用于計(jì)算直接型轉(zhuǎn)換級(jí)聯(lián)型濾波器參數(shù)值的函數(shù)。5 結(jié)束語本次課程設(shè)計(jì)，我的課題是利用Flattopwin窗函數(shù)設(shè)計(jì)FIR濾波器對(duì)音樂信號(hào)濾波去噪。在經(jīng)歷了一個(gè)寒假的放松，腦袋里那點(diǎn)不多的濾波器的只是，早已魂飛天外。但是在實(shí)習(xí)的過程中，通過溫習(xí)課本以及通過頻譜圖，對(duì)濾波器的性能以及功能有了更加深刻的了解。每次課程設(shè)計(jì)都是一次歷練，是對(duì)知識(shí)的查漏補(bǔ)缺，是一次給自己深刻反思的機(jī)會(huì)。每一次實(shí)習(xí)，都可以意識(shí)到自己所欠缺的，意識(shí)到在實(shí)踐面前理論是如何的單薄。知識(shí)的創(chuàng)造永遠(yuǎn)是為了人類的發(fā)展

20、而服務(wù)的，只存在于理論的知識(shí)是無法創(chuàng)造它所應(yīng)有的價(jià)值，而進(jìn)步卻又來源于實(shí)踐和人民的需求，我希望我們用時(shí)間去學(xué)習(xí)的東西能發(fā)揮它應(yīng)有的價(jià)值，甚至于超越。每一次實(shí)踐都可以學(xué)到許多課本上所沒有的只是，也了解到以前所沒有注意到的方面，而在時(shí)間時(shí)，不斷鍛煉的對(duì)立思考的能力更是能讓我們終生受益。相對(duì)與知識(shí)，意識(shí)總是顯得尤為重要。而在實(shí)踐中，重溫的不僅僅是知識(shí)更是，我們?cè)谟谝院笫堑囊庾R(shí)，去獨(dú)立的思考問題，去查閱書籍，去探尋奧妙。而這些是我們平時(shí)所不能學(xué)到的。我總是希望自己能從理論走向時(shí)間，從書本走向課本，希望自己能用自己所學(xué)的為這個(gè)發(fā)展的社會(huì)的建設(shè)獻(xiàn)上自己的一點(diǎn)微薄之力。參考文獻(xiàn)1張志涌精通MATLAB M6

21、.5版北京：北京航空航天大學(xué)出版社，20032約翰··普羅克斯數(shù)字信號(hào)處理M西安：西安交通大學(xué)出版社，20093張小虹信號(hào)系統(tǒng)與數(shù)字信號(hào)處理M第1版西安：西安電子科技出版社，20024謝德芳數(shù)字信號(hào)處理M北京：科學(xué)出版社，20055郝建軍數(shù)字通信M第2版北京：北京郵電大學(xué)出版社，20106張威MATLAB 基礎(chǔ)與編程入門M西安：西安電子科技大學(xué)出版社，2010附錄一： 音樂信號(hào)濾波去噪使用FLARTOPWIN濾波器%程序功能：在Matlab中，用窗口設(shè)計(jì)法設(shè)計(jì)FIR濾波器%程序作者：蔣敏%最后修改日期：2016-3-1程序一-(在音樂信號(hào)中加入噪聲干擾)-x,fs,bits

22、=wavread('e:yinyue.wav'); % 輸入?yún)?shù)為文件的全路徑和文件名，輸出的第一個(gè)參數(shù)是每個(gè)樣本的值，fs是生成該波形文件時(shí)的采樣率，bits是波形文件每樣本的編碼位數(shù)。sound(x,fs,bits); % 按指定的采樣率和每樣本編碼位數(shù)回放N=length(x); % 計(jì)算信號(hào)x的長度fn=2200; % 單頻噪聲頻率t=0:1/fs:(N-1)/fs; % 計(jì)算時(shí)間范圍，樣本數(shù)除以采樣頻率x=x' y=x+0.1*sin(fn*2*pi*t); %在原信號(hào)中加入噪聲sound(y,fs,bits); % 應(yīng)該可以明顯聽出有尖銳的單頻嘯叫聲程序二-

23、(繪制加噪前后音樂信號(hào)的時(shí)間波形和頻譜)- X=abs(fft(x); Y=abs(fft(y); % 對(duì)原始信號(hào)和加噪信號(hào)進(jìn)行fft變換,取幅度譜X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分deltaf=fs/N; % 計(jì)算頻譜的譜線間隔f=0:deltaf:fs/2-deltaf; % 計(jì)算頻譜頻率范圍figure(); %繪制加噪前后音樂信號(hào)的時(shí)間波形和頻譜subplot(221);plot(t,x);title('原始音樂信號(hào)');xlabel('時(shí)間（單位 ： s ）');ylabel('幅度'); grid on;

24、subplot(223);plot(t,y);title('加入干擾信號(hào)的音樂信號(hào)');xlabel('時(shí)間（單位 ： s ）');ylabel('幅度');grid on;subplot(222);plot(f,X);title('音樂信號(hào)幅度譜');xlabel('頻率（單位 ： Hz ）');ylabel('幅度譜');grid on;subplot(224);plot(f,Y);axis(0,4000,0,400); title('加入干擾信號(hào)的音樂信號(hào)幅度譜');xlabe

25、l('頻率（單位 ： Hz ）');ylabel('幅度譜');grid on;程序三-(設(shè)計(jì)濾波器)-fpd=1825;fsd=2075;fsu=2325;fpu=2575;Rp=1;As=20; % 帶阻濾波器設(shè)計(jì)指標(biāo)fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu); % 計(jì)算上下邊帶中心頻率，和頻率間隔wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; % 將Hz為單位的模擬頻率換算為rad為單位的數(shù)字頻率wsd=fsd/fs*2*pi;wsu=fs

26、u/fs*2*pi;wpd=fpd/fs*2*pi; wpu=fpu/fs*2*pi; M=ceil(8.5*pi/dw)+1 ; % 計(jì)算窗設(shè)計(jì)該濾波器時(shí)需要的階數(shù)n=0:M-1; % 定義時(shí)間范圍w_par=(flattopwin (M); % 產(chǎn)生M階的Flattopwin窗hd_bs=IDEAL_LP(wcd,M)+IDEAL_LP(pi,M)-IDEAL_LP(wcu,M); % 調(diào)用自編函數(shù)計(jì)算理想帶阻濾波器的脈沖響應(yīng)h_bs=w_par'.*hd_bs; % 用窗口法計(jì)算實(shí)際濾波器脈沖響應(yīng)db,mag,pha,grd,w=FREQZ_M(h_bs,1); % 調(diào)用自編函數(shù)

27、計(jì)算濾波器的頻率特性dw=2*pi/1000;RP=-min(db(1:wpd/dw+1,ceil(wpu/dw)+1:end) %計(jì)算濾波器的指標(biāo)檢查是否達(dá)標(biāo)AS=-max(db(wsd/dw+1:wsu/dw+1)figure() ;subplot(2,2,1);plot(w/pi,db);title('濾波器幅度響應(yīng)圖');xlabel('w/pi');ylabel('db');axis(0,1,-70,10);line(0,1,-As,-As,'color','r','linestyle',

28、'-','LineWidth',2);line(0,1,-Rp,-Rp,'color','r','linestyle','-','LineWidth',2);line(wsu/pi,wsu/pi,-70,10,'color','r','linestyle','-','LineWidth',2);line(wsd/pi,wsd/pi,-70,10,'color','r',&#

29、39;linestyle','-','LineWidth',2);subplot(2,2,2);plot(w/pi,mag);title('濾波器幅度響應(yīng)圖');xlabel('w/pi');ylabel('幅度mag'); grid on;subplot(2,2,3);plot(w/pi,pha);title('濾波器相位響應(yīng)圖');xlabel('w/pi');ylabel('相位pha'); grid on;subplot(2,2,4);stem(n,h

30、_bs);title('濾波器脈沖響應(yīng)圖');xlabel('n');ylabel('h(n)');grid on;程序四-(濾波)-y_fil=filter(h_bs,1,y); % 用設(shè)計(jì)好的濾波器對(duì)y進(jìn)行濾波Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2); % 計(jì)算頻譜取前一半figure();subplot(321);plot(t,x);title('原始音樂信號(hào)');xlabel('時(shí)間（單位 ： s ）');ylabel('幅度'); grid on;su

31、bplot(323);plot(t,y);title('加入干擾信號(hào)的音樂信號(hào)');xlabel('時(shí)間（單位 ： s ）');ylabel('幅度');grid on;subplot(322);plot(f,X);title('音樂信號(hào)幅度譜');xlabel('頻率（單位 ： Hz ）');ylabel('幅度譜');grid on;subplot(324);plot(f,Y);axis(0,4000,0,400); title('加入干擾信號(hào)的音樂信號(hào)幅度譜');xlabel(

32、'頻率（單位 ： Hz ）');ylabel('幅度譜');grid on;subplot(325);plot(t,y_fil); title('濾波后的音樂信號(hào)');xlabel('時(shí)間（單位 ： s ）');ylabel('幅度');grid on;subplot(326);plot(f,Y_fil);title('濾波后的音樂信號(hào)幅度譜');xlabel('頻率（單位 ： Hz ）');ylabel('幅度譜'); grid on;sound (y_fil,fs

33、,bits);C,B,A=DIR2CAS(h_bs,1); %計(jì)算級(jí)聯(lián)型濾波器參數(shù)附錄二：函數(shù)FREQZ_M.M定義：function db,mag,pha,grd,w = freqz_m(b,a);% freqz 子程序的改進(jìn)版本% -% db,mag,pha,grd,w = freqz_m(b,a);% db = 0 到pi弧度區(qū)間內(nèi)的相對(duì)振幅(db)% mag = 0 到pi弧度區(qū)間內(nèi)的絕對(duì)振幅% pha = 0 到pi弧度區(qū)間內(nèi)的相位響應(yīng)% grd = 0 到pi弧度區(qū)間內(nèi)的群遲延% w = 0 到pi弧度區(qū)間內(nèi)的501個(gè)頻率樣本向量% b = Ha(z)的分子多項(xiàng)式系數(shù)(對(duì)FIR b=h)% a = Ha(z)的分母多項(xiàng)式系數(shù)(對(duì) FIR: a=1)%H,w = freqz(b,a,1000,'whole'); H = (H(1:1:501)