1、 李柳 音樂信號濾波去噪使用脈沖響應不變法設計的切比雪夫I型濾波器 第頁共19頁 音樂信號濾波去噪使用脈沖響應不變法設計的切比雪夫I型濾波器學生姓名：李柳 指導老師：黃紅兵摘 要 本課程設計主要內容是使用脈沖響應不變法設計的切比雪夫I型濾波器，對一段音樂信號進行濾波去噪處理并根據濾波前后的波形和頻譜分析濾波性能。本課程設計仿真平臺為MATLAB7.0，開發工具是M語言編程。首先下載一段音樂，并人為加入一單頻噪聲，然后對信號進行頻譜分析以確定所加噪聲頻率，并設計濾波器進行濾波去噪處理，最后比較濾波前后的波形和頻譜并進行分析。由分析結果可知，濾波器后的音樂信號與原始信號基本一致，即設計的濾波器能夠

2、去除信號中所加單頻噪聲，達到了設計目的。 關鍵詞 濾波去噪；脈沖響應不變法；切比雪夫I型濾波器；MATLAB 1 引 言此次課程設計主要是在網上采集一段8000Hz，8位的單聲道PCM格式音樂信號，并繪制波形觀察其時域和頻域的波形圖，再在MATLAB平臺上，將該音樂信號進行濾波去噪處理，對比濾波前和濾波后的時域和頻域的波形圖，根據結果和學過的理論得出合理的結論。1.1 課程設計目的課程設計有利于我們對基礎知識的理解，并將所學的知識應用起來，此次課程設計用到Matlab，數字信號處理，以及辦公軟件Visio等知識，平時總是在分析濾波器，其實并不太理解濾波器跟我們的生活有什么聯系，而課程設計要求我

3、們自己動手操作，從原始信號的采集到加入噪聲之后的信號到使用我們自己設計的濾波器對加噪信號進行濾波處理之后恢復出原始信號，這個過程讓我們真正理解我們所學的知識在我們生活中的用處，從而讓增強我們對這門學科以及我們專業的興趣。另外此次課程設計也有利于邏輯思維的鍛煉，數字信號處理課程設計是在學生完成數字信號處理和MATLAB的結合后的基本實驗以后開設的。本課程設計的目的是為了讓學生綜合數字信號處理和MATLAB并實現一個較為完整的小型濾波系統。這一點與驗證性的基本實驗有本質性的區別。平時的學習都是分模塊進行，并沒有系統的自己一個人獨立完成設計到操作的過程，這樣系統的設計正好鍛煉了我們這方面的能力。開設

4、課程設計課程的主要目的是通過系統設計、軟件仿真、程序安排與調試、寫課程設計報告等步驟，使學生初步掌握工程設計的具體步驟和方法，提高分析問題和解決問題的能力，提高實際應用水平。 1.2課程設計的要求（1）濾波器指標必須符合工程實際。（2）設計完后應檢查其頻率響應曲線是否滿足指標。（3）處理結果和分析結論應該一致，而且應符合理論。（4）獨立完成課程設計并按要求編寫課程設計報告書。1.3設計平臺此次課程設計在MATLAB平臺下進行設計，MATLAB（矩陣實驗室）是MATrix LABoratory的縮寫，是一款由美國The MathWorks公司出品的商業數學軟件。MATLAB是一種用于算法開發、數

5、據可視化、數據分析以及數值計算的高級技術計算語言和交互式環境。除了矩陣運算、繪制函數/數據圖像等常用功能外，MATLAB還可以用來創建用戶界面及與調用其它語言（包括C，C+和FORTRAN）編寫的程序。MATLAB可以進行矩陣運算、繪制函數和數據、實現算法、創建用戶界面、連接其他編程語言的程序等，主要應用于工程計算、控制設計、信號處理與通訊、圖像處理、信號檢測、金融建模設計與分析等領域。高效的數值計算及符號計算功能，能使用戶從繁雜的數學運算分析中解脫出來，且具有完備的圖形處理功能，實現計算結果和編程的可視化，友好的用戶界面及接近數學表達式的自然化語言，使學者易于學習和掌握，功能豐富的應用工具箱

6、(如信號處理工具箱、通信工具箱等) ，為用戶提供了大量方便實用的處理工具。2 設計原理在網上采集一段音樂信號，繪制波形并觀察其頻譜，給定相應技術指標，用脈沖響應不變法設計的一個滿足指標的切比雪夫I型IIR濾波器，并對該信號進行濾波去噪處理，之后比較濾波前后的波形和頻譜并進行分析。2.1 IIR濾波器IIR濾波器具有無限長脈沖響應，因此能夠與模擬濾波器相匹敵；一般來說，所有的模擬濾波器都有無限長脈沖響應。因此，IIR濾波器設計的基本方法是利用復值映射將大家熟知的模擬濾波器變換為數字濾波器。這個基本方法稱為A/D（模擬-數字）濾波器變換。IIR數字濾波器采用遞歸型結構，即結構上帶有反饋環路。IIR

7、濾波器運算結構通常由延時、乘以系數和相加等基本運算組成，可以組合成直接型、級聯型、并聯型三種結構形式，都具有反饋回路。通常有IIR數字濾波器的直接和間接設計法，所謂模擬濾波器設計數字濾波器，第一步將設計的歸一化樣本模擬低通濾波器經模擬模擬頻帶變換法轉換成模擬低通濾波器，第二步，然后數字化采用脈沖響應不變法成各相應頻帶的數字濾波器2。2.2切比雪夫I型濾波器 切比雪夫濾波器，又名“車比雪夫濾波器”，是在通帶或阻帶上頻率響應幅度等波紋波動的濾波器。切比雪夫濾波器來自切比雪夫分布，以“切比雪夫”命名，是用以紀念俄羅斯數學家巴夫尼提列波維其切比雪夫。切比雪夫濾波器在過渡帶比巴特沃斯濾波器的衰減快，但頻

8、率響應的幅頻特性不如后者平坦。切比雪夫濾波器和理想濾波器的頻率響應曲線之間的誤差最小，但是在通頻帶內存在幅度波動。I型切比雪夫濾波器在通帶（或稱“通頻帶”）上頻率響應幅度等波紋波動的濾波器稱為“I型切比雪夫濾波器；n階第一類切比雪夫濾波器的幅度與頻率的關系可用下列公式表示： （2.1）其中： 而 是濾波器在截止頻率：的放大率 是 n階切比雪夫多項式。 （2.2） （2.3）其中 或: （2.4） （2.5）切比雪夫濾波器的階數等于此濾波器的電子線路內的電抗元件數。切比雪夫濾波器的幅度波動 = 分貝當 ,切比雪夫濾波器的幅度波動= 3分貝。如果需要幅度在在阻頻帶邊上衰減得更陡峭，可允許在復平面的

9、jw軸上存在零點。但結果會使通頻帶內振幅波動較大，而在阻頻帶內對信號抑制較弱。切比雪夫濾波器存在兩種可能的幅度平方響應的形狀，當N分別為奇數和偶數時如圖2-1所示： (a) N為奇數 （b）N為偶數 圖21 切比雪夫I型濾波器的幅度特性 2.3脈沖響應不變法脈沖響應不變法是從濾波器的脈沖響應出發，使數字濾波器的單位脈沖響應序列h(n) 逼近模擬濾波器的沖擊響應ha(t), 使h(n)正好等于ha(t)的抽樣值，即滿足 （2.6） 其中，T為采樣周期2。如以Ha(s)及H（z）分別表示ha(t)的拉氏變換及h(n)的z變換，設 且設只有單級階點，且假定的分母階次NM，則根據采樣序列z變換與模擬信

10、號拉氏變換的關系，得: （2.7） 上式表明，采用脈沖響應不變法將模擬濾波器變換為數字濾波器時，它所完成的S平面到Z平面的變換，正是以前討論的拉氏變換到Z變換的標準變換關系，即首先對Ha(s)作周期延拓，然后再經過z=est的映射關系映射到Z平面上2。 z=est的映射關系表明，S平面上每一條寬為2/T的橫帶部分，都將重疊地映射到Z平面的整個全部平面上。每一橫帶的左半部分映射到Z平面單位圓以內，每一橫帶的右半部分映射到Z平面單位圓以外，j軸映射在單位圓上，但j軸上的每一段2/T都對應于繞單位圓一周，如下圖所示：(a)s平面 (b) z平面圖2-2 脈沖響應不變法的映射關系應當指出，Z=est的

11、映射關系反映的是Ha(s)的周期延拓與H（z）的關系，而不是Ha(s)本身與H（z）的關系，因此，使用脈沖響應不變法時，從Ha(s)到H(z)并沒有一個由S平面到Z平面的簡單代數映射關系，即沒有一個s=f(z)的代數關系式。另外，數字濾波器的頻響也不是簡單的重現模擬濾波器的頻響，而是模擬濾波器頻響的周期延拓，周期為S=2/T=2fs，即 （2.8）脈沖響應不變法特別適用于用部分分式表達的傳遞函數，模擬濾波器的傳遞函數若只有單階極點，且分母的階數高于分子階數NM，則可表達為部分分式形式：其拉氏反變換為： 其中u(t)為單位階躍函數。對ha(t)采樣就得到數字濾波器的單位脈沖響應序列，再對h(n)

12、取Z變換，得到數字濾波器的傳遞函數： （2.9）第二個求和為等比級數之和，要收斂的話，必有，所以有: （2.10）比較部分分式形式的Ha(s)和上式H(z)可以看到，S平面上的極點s=si,變換到Z平面上是極點,而Ha(s)與H(z)中部分分式所對應的系數不變。如果模擬濾波器是穩定的，則所有極點si都在S左半平面，即Resi0,那么變換后H(z)的極點也都在單位圓以內，即，因此數字濾波器保持穩定。雖然脈沖響應不變法能保證S平面與Z平面的極點位置有一一對應的代數關系，但這并不是說整個S平面與Z平面就存在這種一一對應的關系，特別是數字濾波器的零點位,與S平面上的零點就沒有一一對應關系，而是隨著的極

13、點與系數的不同而不同。 是的周期延拓（周期為fs），因并不是帶限，即在超過頻率部分并不為0，所以就產生了混迭。當為低通或帶通濾波器時， 越大，則的下一周期相隔越遠，混迭也就越小。當為帶阻或高通濾波器時，在超過頻率部分全為通帶，這樣就不滿足抽樣定理，發生了完全的混迭，所以脈沖響應不變法不能設計帶阻或高通濾波器。3設計步驟3.1 設計流程圖認真閱讀老師下發的任務書，上網搜集查找相關材料，大致的畫出設計流程圖，即音樂信號濾波去噪使用脈沖不變響應法設計的切比雪夫I型濾波器的設計流程如圖3-1所示： 圖3-1 設計流程圖3.2 錄制音樂信號 在網上下載一段音樂信號，利用windows下的錄音機工具通過點

14、擊錄音機界面上的文件屬性，在彈出的對話框中選擇馬上轉換，設置錄音格式為PCM編碼，屬性為8kHz，單聲道，這樣的設置在相同的時間下數據最少，處理時間也就最短。點擊文件，保存，存為全英文的文件名，后綴為“.wav”。如“liliu.wav”。并將其放入MATLAB安裝盤下的work文件夾中就好了。如下圖所示： 圖3-2 音樂信號格式轉換圖 8000Hz的采樣率是指：聲音信號在“模數”轉換過程中單位時間內采樣8000次。采樣值是指每一次采樣周期內聲音模擬信號的積分值。對于單聲道聲音文件，采樣數據為八位的短整數，wave文件是計算機領域最常用的數字化聲音文件格式之一，它是微軟專門為Windows系統

15、定義的波形文件格式（Waveform Audio），由于其擴展名為*.wav。PCM也被稱為 脈碼編碼調制。PCM中的聲音數據沒有被壓縮，如果是單聲道的文件，采樣數據按時間的先后順序依次存入,它的基本組織單位是BYTE(8bit)。 采集信號完成以后，利用MATLAB畫出原始信號在時域和頻域的波形圖，編寫M文件如下： x,fs,bits=wavread(liliu.wav); % 輸入參數為文件的全路徑和文件名，輸出的第一個參數是每個樣本的值，fs是生成該波形文件時的采樣率，bits是波形文件每樣本的編碼位數。sound(x,fs,bits); % 按指定的采樣率和每樣本編碼位數回放N=len

16、gth(x); % 計算信號x的長度t=0:1/fs:(N-1)/fs; % 計算時間范圍，樣本數除以采樣頻率X=abs(fft(x); % 對原始信號進行fft變換,取幅度譜X=X(1:N/2); % 截取前半部分deltaf=fs/N; % 計算頻譜的譜線間隔f=0:deltaf:fs/2-deltaf; % 計算頻譜頻率范圍得到的原始信號的時域與頻域的波形圖如下圖所示： 圖3-3 原始音樂信號時域與頻域的波形畫出音樂信號的時域波形；然后對音樂號進行快速傅里葉變換，得到信號的頻譜特性，從原始音樂信號的幅度譜中可以看出，音樂信號的頻率主要集中在01000Hz ，2000Hz以后的頻率段幾乎沒

17、有能量分布。采集到合適格式的音樂信號之后，對原始音樂信號進行加噪處理，然后對加入單頻率的噪聲的音樂信號進行fft變換,取幅度譜，畫出加噪的音樂信號在時域和頻域的波形，并對比進行分析，主要程序如下：x=x; y=x+0.1*sin(fn*2*pi*t); %加噪信號sound(y,fs,bits); %加噪信號聲音Y=abs(fft(y); %對原始信號和加噪信號進行fft變換,取幅度譜 Y=Y(1:N/2); % 截取前半部分deltaf=fs/N; % 計算頻譜的譜線間隔f=0:deltaf:fs/2-deltaf;% 計算頻譜頻率范圍運行程序的得到的加噪信號的時域與頻域的波形圖如下所示：

18、圖3-4 加噪信號的時域與頻域的波形圖 由于此次課程設計加入的噪聲為fn=2000Hz的單頻信號，運行程序的時候可以明顯的聽到加噪之后的音樂信號對比原始信號有明顯的單頻呼嘯聲，由運行結果圖3-4右下角的頻譜圖也可明顯的看到在頻率為2000Hz時存在單頻噪聲，時域圖也有明顯的變化。3.3 濾波器設計我們這次設計的主要任務是對單頻干擾進行濾波處理，因此設計一個低通濾波器，設計指標如下：取fp=fn-200=1800 fc=fn-50=1950，通帶波紋，阻帶衰減，用脈沖響應不變法設計的切比雪夫I型濾波器對單頻干擾進行濾波處理。設計濾波器的重要程序如下：fp=fn-200;fc=fn-50; %定義

19、通帶和阻帶截止頻率Rp=1;As=40; % 定義通帶波紋和阻帶衰減wp=fp/fs*2*pi;ws=fc/fs*2*pi; %計算對應的數字頻率T=1; %定義采樣間隔Omegap=wp/T; Omegas=ws/T; %模數指標變換N,Omegac=cheb1ord(Omegap,Omegas,Rp,As,s) %計算切比雪夫濾波器階數和截止頻率cs,ds=cheby1(N,Rp,Omegac,s);% 計算切比雪夫濾波器系統函數分子分母系數% 第三步 模數濾波器變換 b,a=impinvar(cs,ds,1/T) % 脈沖響應不變法得到數字濾波器系統函數分子分母系數db,mag,pha,

20、grd,w=freqz_m(b,a); %驗證濾波器是否達到指定性能delta=1,zeros(1,99);ha=filter(b,a,delta);% 計算脈沖響應在MATLAB中運行濾波器的程序，得到濾波器的特性圖，如下圖所示：圖3-5 濾波器特性圖 由程序運行結果可知wp/pi=0.450 ，ws/pi =0.4875，N =15， Omegac =1.4137,對比圖3-4中濾波器的幅度特性圖可知當w/pi0.4875時，dB都大于40，所以利用脈沖響應不變法設計的切比雪夫濾波器符合要求。3.4 信號濾波處理在3.3節中已經用使用脈沖響應不變法設計的切比雪夫I型濾波器，現在只需用自己設

21、計的切比雪夫I型濾波器對采集的信號進行濾波，在Matlab中， IIR濾波器利用函數filter對含噪信號y進行濾波去噪處理，得到濾波后的信號存儲為y1，并計算y1的頻譜。部分程序如下：y1=filter(b,a，y); % 用設計好的濾波器對y進行濾波處理,得到時域圖Y1=fft(y1);Y1=Y1(1:length(Y1)/2); % 計算頻譜取前一半 對比濾波前后的時域與頻域的波形圖3-5，具體分析濾波是否達到要求。 圖3-6 濾波前后的時域與頻域的波形對比圖我們觀察到圖3-5濾波前后音樂信號的波形對比圖，發現在時域波形中濾波后的波形與加干擾噪聲后有明顯的變化，與原始音樂信號幾乎相同，說

22、明設計的濾波器達到要求。而且在頻域波形中，很明顯地反應出設計的濾波器濾去了我們采集的音樂信號中的噪聲，即濾去了頻率高于2000Hz的信號。3.5 結果分析綜合上述的實驗步驟及其實驗得出的結果圖，在圖34中我們可以發現加入的2000Hz單頻噪聲是直接累加在原始音樂信號上的，屬于加性噪聲，在圖3-5 濾波器特性圖中可以得出，我們利用脈沖響應不變法設計的切比雪夫I濾波器屬于低通濾波器，能將2000Hz以上的噪聲濾掉，分析圖36濾波前后的時域與頻域的波形對比圖，我們就能明顯看出我們利用脈沖響應不變法設計的切比雪夫I濾波器濾去了我們采集的音樂信號中的噪聲，達到了設計要求。3.6 濾波器結構設計 IIR濾

23、波器可以組合成直接型、級聯型、并聯型三種基本結構形式，都具有反饋回路。根據設計的濾波器的原始代碼的運行結果，我們可以得知濾波器的階數N15，相對來說濾波器的階數不是特別大，這里我們采用直接型的結構。濾波器的結構圖如下圖所示：圖36 濾波器結構圖由濾波器的源代碼的運行結果可以知道系統函數分子分母系數b和a的值，如下圖所示： 圖37數字濾波器系統函數分子分母系數根據系統函數的分子分母系數以及IIR濾波器直接型結構的特點很容易可以畫圖36濾波器結構圖，N階濾波器需要N個延時單元，左側從上到下是分母系數的負數,右側是分子系數。4出現的問題及解決方法這次的課程設計中我們設計的濾波器對音樂信號的濾噪處理或

24、多或少會因為我們不仔細和實踐能力不足會遇到各種問題，當遇到問題時我們運用自己學到的知識以及同學和老師的幫助將其解決，不斷完善和修改。設計中出現的問題及解決方法有以下幾個方面： （1）在設計的過程中，剛開始從網上下載的音樂信號是立體聲，不符合課程設計的要求，經過和同學的討論，利用WINDOWS下的錄音機工具通過點擊錄音機界面上的文件屬性，在彈出的對話框中選擇馬上轉換，設置錄音格式為PCM編碼，屬性為8kHz，單聲道這樣就能符合原始信號的要求了。 （2）加噪后的信號的幅度譜圖，單頻噪聲幅度遠遠高于其他頻譜分量，導致其他頻譜分量幅度過低而看不清楚加噪前后的頻譜關系，通過使用axis調整縱坐標的范圍即

25、可解決。（3）其次在設濾波器的參數時取fp=fn-200=1800 fc=fn-50=1950，通帶波紋Rp=1dB阻帶衰減As=16dB，濾波后發現不能完全將噪聲濾掉，仍然存在噪聲，在設置的通帶截止頻率和阻帶截止頻率的差值太大就是太小，都不能完全將噪聲濾掉，經過一次次調整與修改最終將參數設置為fp=fn-200=1800 fc=fn-50=1950，通帶波紋Rp=1dB，阻帶衰減As=40dB，這樣就能更好的將噪聲濾除，基本達到指標要求。（4）剛開始設計好濾波器是，不知道怎樣才是符合要求，將設計好的濾波器拿給老師看，老師說在Wp/pi,Ws/pi以及Rp,As，出各畫一條直線能更好的看出濾波

26、器是否滿足要求，即，當w/pi大與ws/pi時都大于As=40dB就說明滿足要求了。（5）最后就是在使用Visio軟件時的問題，由于是第一次使用，比較陌生，最開始使用的線是尾端帶箭頭的，這樣每次畫一段就需要兩次修改，最后仔細的在工具箱中尋找，發現并使用了中間帶箭頭的線，這樣畫起圖來方便很多。5 結束語通過這段時間的親身經歷，我感覺自己學到了：收集、整理資料、共同協作、分析及處理問題等許多方面的知識。在做設計的過程中，我有許多不懂得地方，在老師的指導下我一步步的解決問題完成論文，在完成過程中老師指導我去怎么選擇資料，如何去利用網絡資源，在這個學習的過程中，我了解到MATLAB的實用價值，同時利用

27、MATLAB簡單編程的特點，可以實現較為復雜的系統，因此MATLAB在通信系統仿真方面具有強大的功能和優越性，在課程設計的過程中明顯覺得自己對MATLAB軟件的不熟練，很多函數不會使用，而這次的課程設計讓我深刻認識到MATLAB的實用性，同時也更加熟練的使用MATLAB，掌握了MATLAB的更多用法。同時這次課程設計是數字信號處理的，在設計的過程中明顯發現對課本知識的不熟練，用到的很多課本知識仍要去翻閱課本，這讓我認識到在以后的學習中，應該更加認真更加深入，不能淺嘗輒止。在以后的學習或者其他設計中，應該提前做好充足的準備，對要做的東西有足夠的了解之后再開始進行，其次以后在工作中要腳踏實地，在學

28、術上要嚴謹，在思維上要活躍，在學業上要勤奮刻苦，不能懵懂的就以為自己已經會了。最后我真誠感謝這期間老師給予我的全力幫助，細心指導以及對我的嚴格要求，是老師在我遇到問題時，不辭辛苦幫我解決，感謝老師在設計和任務安排上長時間的指導。參考文獻1 張圣勤.MATLAB7.0實用教程M.北京：機械工程出版社,2006.P140-180 2 程佩青.數字信號處理教程M.北京：清華大學出版社,2014,P269-4633 張志美. MATLAB完全自學手冊M.電子工業出版社,2013,P142-405 4 邵毅,張明玉,溫艷,基于MATLAB的IIR數字濾波器的設計J.宿州學院學報,2009年02期 5 王

29、艷文,楊楠, 基于MATLAB/脈沖響應不變法的切比雪夫I型數字濾波器設計J. 新鄉醫學院生命科學技術學院.2013年18期 附錄：音樂信號濾波去噪設計源程序清單% 程序名稱：liliu.m% 程序功能：音樂信號濾波去噪使用脈沖響應不變法設計的切比雪夫I型濾波器% 程序作者：李柳% 最后修改日期：2015-3-13x,fs,bits=wavread(liliu.wav); % 輸入參數為文件的全路徑和文件名，輸出的第一個參數是每個樣本的值，fs是生成該波形文件時的采樣率，bits是波形文件每樣本的編碼位數。sound(x,fs,bits); % 按指定的采樣率和每樣本編碼位數回放N=lengt

30、h(x); % 計算信號x的長度fn=2000; % 單頻噪聲頻率，此參數可改t=0:1/fs:(N-1)/fs; % 計算時間范圍，樣本數除以采樣頻率x=x; y=x+0.1*sin(fn*2*pi*t);sound(y,fs,bits); %回放加噪音樂X=abs(fft(x); Y=abs(fft(y); % 對原始信號和加噪信號進行fft變換,取幅度譜X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分deltaf=fs/N; % 計算頻譜的譜線間隔f=0:deltaf:fs/2-deltaf;% 計算頻譜頻率范圍fn=2000; % 單頻噪聲頻率fs=8000;fp=fn

31、-200;fc=fn-50; %定義通帶和阻帶截止頻率Rp=1;As=40; % 定義通帶波紋和阻帶衰減wp=fp/fs*2*pi;ws=fc/fs*2*pi; %計算對應的數字頻率T=1; %定義采樣間隔% 第一步 模數指標變換Omegap=wp/T; Omegas=ws/T; %截止頻率線性變換% 第二步 根據指標設計模擬濾波器N,Omegac=cheb1ord(Omegap,Omegas,Rp,As,s) %計算濾波器階數和截止頻率cs,ds=cheby1(N,Rp,Omegac,s);% 計算濾波器系統函數分子分母系數% 第三步 模數濾波器變換b,a=impinvar(cs,ds,1/

32、T) % 脈沖響應不變法得到數字濾波器系統函數分子分母系數db,mag,pha,grd,w=freqz_m(b,a); %驗證濾波器是否達到指定性能delta=1,zeros(1,99);ha=filter(b,a,delta);% 計算脈沖響應Rp = -min(db(1:fix(wp/(2*pi/1000)+1) % 實際的通帶波動As = -round(max(db(fix(ws/(2*pi/1000)+1:501) % 最小阻帶衰減y1=filter(b,a,y) % IIR濾波器對信號signal進行濾波處理，其中b,a為上面設計好的濾波器參數。Y1=abs(fft(y1);Y1=Y

33、1(1:length(Y1)/2); % 計算頻譜取前一半figure(1)subplot(221);plot(t,x) % 繪制截取后的語音信號title(原始音樂信號);xlabel(時間(單位：s);ylabel(幅度);axis tightsubplot(223);plot(t,y) % 繪制截取后的語音信號title(加入單頻干擾后的音樂信號);xlabel(時間(單位：s);ylabel(幅度);axis tightsubplot(222);plot(f,abs(X) % 繪制頻譜圖title(原始音樂信號幅度譜圖);xlabel(頻率（單位：Hz）);ylabel(幅度譜) ;a

34、xis tightsubplot(224);plot(f,abs(Y) % 繪制頻譜圖title(加入干擾后的音樂信號幅度譜圖);xlabel(頻率（單位：Hz）);ylabel(幅度譜) ;axis (0 4000 0 400)figure(2)Subplot(221);plot(w/pi,db);axis(0 1 -50 20);title(濾波器幅度響應);xlabel(w/pi);ylabel(dB);grid on;line(wp/pi,wp/pi,-50,20,Color,r,LineWidth,2,LineStyle,-);line(Omegac/pi,Omegac/pi,-50,20,Color,r,LineWidth,2,LineStyle,-);l