使用脈沖響應(yīng)不變法設(shè)計(jì)切比雪夫I型IIR濾波器

語(yǔ)音信號(hào)濾波去噪—脈沖響應(yīng)不變法設(shè)計(jì)的并聯(lián)型切比雪夫I

型濾波器學(xué)生姓名：指導(dǎo)老師：摘要本課程設(shè)計(jì)主要內(nèi)容是使用脈沖響應(yīng)不變法設(shè)計(jì)的并聯(lián)型切比雪夫I型IIR濾波器，對(duì)一段含噪語(yǔ)音信號(hào)進(jìn)行濾波去噪處理并根據(jù)濾波前后的波形和頻譜分析濾波性能。本課程設(shè)計(jì)仿真平臺(tái)為MATLAB7.0,開(kāi)發(fā)工具是M語(yǔ)言編程。首先利用錄音工具錄制一段語(yǔ)音信號(hào)，并人為加入一單頻噪聲，然后對(duì)信號(hào)進(jìn)行頻譜分析以確定所加噪聲頻率，并設(shè)計(jì)濾波器進(jìn)行濾波去噪處理，最后比較濾波前后的波形和頻譜并進(jìn)行分析。由分析結(jié)果可知，濾波器后的語(yǔ)音信號(hào)與原始信號(hào)基本一致，即設(shè)計(jì)的切比雪夫I型IIR濾波器能夠去除信號(hào)中所加單頻噪聲，達(dá)到了設(shè)計(jì)目的。關(guān)鍵詞濾波去噪；脈沖響應(yīng)不變法；切比雪夫I型；IIR濾波器；MATLAB1引言本課程設(shè)計(jì)主要解決在含噪情況下對(duì)語(yǔ)音信號(hào)的濾波去噪處理，處理時(shí)采用的是利用脈沖響應(yīng)不變法設(shè)計(jì)的切比雪夫I型IIR濾波器。1.1課程設(shè)計(jì)目的（1）熟悉使用MATLAB語(yǔ)言來(lái)解決一些簡(jiǎn)單的課程問(wèn)題。（2）鍛煉動(dòng)手能力及獨(dú)立思考并解決問(wèn)題的能力（3）增加對(duì)書(shū)本知識(shí)的理解以及掌握，如：IIR濾波器的工作原理，切比雪夫?yàn)V波器的設(shè)計(jì)，脈沖響應(yīng)法的應(yīng)用等等；（4）對(duì)錄制好的音頻行加噪和過(guò)濾，比較前后語(yǔ)音變化;

1.2課程設(shè)計(jì)內(nèi)容（1）使用錄音軟件;錄制一段語(yǔ)音信號(hào)，要求為PCM編碼格式單聲道8K8位7k/s的語(yǔ)音信號(hào)，繪制波形并觀察頻譜特點(diǎn)；（2）加上噪聲并繪制加噪信號(hào)波形及觀察其頻譜特點(diǎn)；（3）設(shè)計(jì)濾波器，繪制出濾波器的頻譜并檢驗(yàn)濾波器是否合格；（4）最后使用設(shè)計(jì)好的濾波器對(duì)噪聲信號(hào)進(jìn)項(xiàng)過(guò)濾，繪制過(guò)濾后的波形及觀察濾波前后頻譜變化；（5）繪制濾波器的并聯(lián)型結(jié)構(gòu)圖；2設(shè)計(jì)原理錄制一段語(yǔ)音信號(hào)，繪制波形并觀察其頻譜，給定相應(yīng)技術(shù)指標(biāo)，用脈沖響應(yīng)不變法設(shè)計(jì)的一個(gè)滿足指標(biāo)的切比雪夫I型IIR濾波器，對(duì)該語(yǔ)音信號(hào)進(jìn)行濾波去噪處理，比較濾波前后的波形和頻譜并進(jìn)行分析。2.1IIR濾波器IIR濾波器具有無(wú)限長(zhǎng)脈沖響應(yīng)，一般來(lái)說(shuō)，所有的模擬濾波器都有無(wú)限長(zhǎng)脈沖響應(yīng)。因此，IIR濾波器設(shè)計(jì)的基本方法是利用復(fù)值映射將大家熟知的模擬濾波器變換為數(shù)字濾波器。這種IIR濾波器設(shè)計(jì)的基本方法存在途徑之一［2］如下：期望的IIR期望的IIR濾波器2.2切比雪夫I型濾波器切比雪夫I型濾波器是電子濾波器的一種。它的特點(diǎn)是在通帶具有等波紋響應(yīng)。其振幅平方函數(shù)如2-1式所示：(2-1)IHa(jO)|2=

式中，N為整數(shù)，稱為濾波器的階，已是通帶波紋因子(它與Rp有關(guān))。其中N階切比雪夫多項(xiàng)式如2-2式所示：(2-1)(2-2)cos(Ncos-1(x)),0<x<1

cosh(cosh-1(x)),1<x<x(2-2)其中x=—Qc切比雪夫?yàn)V波器存在兩種可能的幅度平方響應(yīng)的形狀，當(dāng)N分別為基數(shù)和偶數(shù)時(shí)如圖2-1所示：圖2-2切比雪夫I型濾波器的幅度特性2.3脈沖響應(yīng)不變法脈沖響應(yīng)不變法是一種將模擬濾波器轉(zhuǎn)化為數(shù)字濾波器的基本方法。它通過(guò)對(duì)模擬濾波器的單位沖激響應(yīng)h(t)等間隔抽樣來(lái)獲取對(duì)應(yīng)數(shù)字濾波器的單位脈沖響應(yīng)h[k],即h[k]=h(t)lt=kT=h(kT)，其中T是等間隔抽樣。不失一般性地，假定連續(xù)時(shí)間濾波器的系統(tǒng)函數(shù)為：NA(2-3)H(s)=Zkk=1(2-3)(2-4)h(t)=rAkeskt，t2°ck=1I°,t<°對(duì)Tdhc(t)采樣得到的離散時(shí)間濾波器的脈沖響應(yīng)是:(2-4)h[n]=Th(nT)=ZtAesknTdu[n]=ZtA(eskTd)nu[n](25)k=1dcddk=1離散時(shí)間濾波器的系統(tǒng)函數(shù)為:H(z)zNTdAk1一eskTzi-1k=1(2-6)3設(shè)計(jì)步驟語(yǔ)音信號(hào)濾波去噪一一使用脈沖響應(yīng)不變法設(shè)計(jì)的并聯(lián)型切比雪夫IIR濾波器，設(shè)計(jì)流程如下圖：開(kāi)始v,V錄制單聲道8KWAV格式的語(yǔ)音信號(hào)繪制語(yǔ)音信號(hào)的波形圖，分析頻譜特點(diǎn)加入噪聲并繪制波形圖，分析噪聲頻譜變化使用脈沖響應(yīng)不變法設(shè)計(jì)切比雪夫1型IIR濾波器；畫(huà)出濾波器并聯(lián)型結(jié)構(gòu)圖對(duì)噪聲信號(hào)進(jìn)行過(guò)濾，畫(huà)出濾波后頻譜圖，分析濾波前后頻譜變化回放語(yǔ)音信號(hào)結(jié)束vy圖3-1步驟圖3.1錄制語(yǔ)音信號(hào)利用錄音軟件(如下圖3-2)錄制PCM編碼格式單聲道8K8位7k/s的語(yǔ)音信號(hào)，語(yǔ)音內(nèi)容不限;圖3-2錄制語(yǔ)音信號(hào)軟件界面3.2繪制語(yǔ)音信號(hào)波形圖使用MATLAB軟件，將語(yǔ)音信號(hào)放到MATLAB的工作目錄中；使用[x,fs.bits]=wavread(‘ldd.wa函數(shù)讀入語(yǔ)音信號(hào)；源程序如下：[x,fs,bits]=wavread('ldd.wav');%輸入?yún)?shù)為文件的全路徑和文件名，輸出的第一個(gè)參數(shù)是每個(gè)樣本的值，fs是生成該波形文件時(shí)的采樣率，bits是波形文件每樣本的編碼位數(shù)。sound(x,fs,bits);%按指定的采樣率和每樣本編碼位數(shù)回放N=length(x);%計(jì)算信號(hào)x的長(zhǎng)度t=0:1/fs:(N-1)/fs;%計(jì)算時(shí)間范圍，樣本數(shù)除以采樣頻率x=x';X=abs(fft(x));%對(duì)原始信號(hào)和加噪信號(hào)進(jìn)行fft變換，取幅度譜X=X(1:N/2);%截取前半部分

deltaf=fs/N;%計(jì)算頻譜的譜線間隔f=0:deltaf:fs/2-deltaf;%計(jì)算頻譜頻率范圍繪制的波形圖如下:圖3-3原始語(yǔ)音信號(hào)的頻譜圖3.3繪制加噪后信號(hào)波形圖在原始語(yǔ)音信號(hào)上加入單頻干擾信號(hào)，fn=3000Hz，其源程序如下：[x,fs,bits]=wavread('ldd.wav');%輸入?yún)?shù)為文件的全路徑和文件名，輸出的第一個(gè)參數(shù)是每個(gè)樣本的值，fs是生成該波形文件時(shí)的采樣率，bits是波形文件每樣本的編碼位數(shù)。sound(x,fs,bits);%按指定的采樣率和每樣本編碼位數(shù)回放N=length(x);%計(jì)算信號(hào)x的長(zhǎng)度f(wàn)n=2000;%單頻噪聲頻率，此參數(shù)可改t=0:1/fs:(N-1)/fs;%計(jì)算時(shí)間范圍，樣本數(shù)除以采樣頻率x=x';y=x+0.1*sin(fn*2*pi*t);sound(y,fs,bits);%應(yīng)該可以明顯聽(tīng)出有尖銳的單頻嘯叫聲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);sound(y,fs,bits);%應(yīng)該可以明顯聽(tīng)出有尖銳的單頻嘯叫聲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ì)算頻譜頻率范圍通過(guò)頻譜圖可以直觀的看到在3000Hz處有尖銳的單峰，與干擾頻率一直；通過(guò)，通過(guò)頻譜圖可以直觀的看到在3000Hz處有尖銳的單峰，與干擾頻率一直；通過(guò)，sound(y,fs,bits)回放該錄音，可以明顯的聽(tīng)到尖銳的單頻噪鳴聲;Rp=1;As=16;%定義通帶波紋和阻帶衰減3.4Rp=1;As=16;%定義通帶波紋和阻帶衰減將數(shù)字濾波器的設(shè)計(jì)指標(biāo)設(shè)為通帶截止頻率fn=3000,fp=fn-200=2800Hz,阻帶頻率fc=fn-50=2950HZ，通帶波紋Rp=1，阻帶波紋As=16dB；濾波器設(shè)計(jì)源程序如下：%定義通帶和阻帶截止頻率%定義通帶和阻帶截止頻率%計(jì)算對(duì)應(yīng)的數(shù)字頻率wp=(fp/fs)*2*pi;ws=(fc/fs)*2*pi;T=1;%%計(jì)算對(duì)應(yīng)的數(shù)字頻率T=1;OmegaP=wp/T;OmegaS=ws/T;%截止頻率線性變換Ripple=10A(-Rp/20);Attn=10A(-As/20);%轉(zhuǎn)化為絕對(duì)形式[cs,ds]=afd_chb1(OmegaP,OmegaS,Rp,As);%計(jì)算濾波器階數(shù)和截止頻率[b,a]=imp_invr(cs,ds,T);%脈沖響應(yīng)不變法得到濾波器系統(tǒng)函數(shù)分子分母系數(shù)[db,mag,pha,grd,w]=freqz_m(b,a);%驗(yàn)證濾波器是否達(dá)到指定性能delta=[1,zeros(1,99)];ha=filter(b,a,delta);%計(jì)算脈沖響應(yīng)運(yùn)行程序得到的濾波器響應(yīng)圖如下:圖3-5濾波器響應(yīng)圖在命令行分別輸入wp/pi,ws/pi驗(yàn)證濾波器是否滿足設(shè)計(jì)標(biāo)準(zhǔn):>>wp/pians=0.7000>>ws/pians=0.7375

3.5對(duì)噪聲信號(hào)進(jìn)行濾波，觀察前后頻譜變化使用設(shè)計(jì)好的濾波器對(duì)噪聲信號(hào)進(jìn)行過(guò)濾，其源程序如下:y_fil=filter(b,a,y);%IIR濾波器對(duì)信號(hào)y進(jìn)行濾波處理，其中b,a為上面設(shè)計(jì)好的濾波器參數(shù)。y_fil=filter(b,a,y);Y_fil=abs(fft(y_fil));%對(duì)fftY_fil=abs(fft(y_fil));%對(duì)fft變換，取幅度譜Y_fil=Y_fil(1:N/2);%截取前半部分sound(y_fil,fs,bits);程序運(yùn)行結(jié)果如下圖:圖3-6圖3-6語(yǔ)音信號(hào)頻譜分析圖由圖3-6可以看出頻率＞3000Hz被過(guò)濾掉，3000Hz處的單頻尖峰噪聲幅度明顯減弱，使用sound(y_fil,fs,bits)函數(shù)播放過(guò)濾后的語(yǔ)音，可以聽(tīng)出過(guò)濾前后的語(yǔ)音變化，說(shuō)明濾波器設(shè)計(jì)基本符合性能要求。3.6繪制濾波器的并聯(lián)型結(jié)構(gòu)圖設(shè)計(jì)完成濾波器后在窗體命令行輸入[C,B,A]=dir2par(b,a)可以得出并聯(lián)型結(jié)構(gòu)的系數(shù)CBA，如下圖所示：***切比雪夫T濾波器階次=10|>>[CjBjA]=dir2par(bja)C=[]B二TOC\o"1-5"\h\z0.09990.0803-0.2599-0.23240.31130.3352-0.2408-0.30660.08950.1260A二1.00001.11150.90611.00000.68880.75121.00000.00020.64041.0000-0.80520.57031.0000-1.37750.5366圖3-7并聯(lián)型結(jié)構(gòu)系數(shù)圖

根據(jù)并聯(lián)型結(jié)構(gòu)系數(shù)以及濾波器階數(shù)可以畫(huà)出濾波器的結(jié)構(gòu)圖如下:4異常處理（1）錄制語(yǔ)音信號(hào)的時(shí)候不了解音頻的編碼方式的區(qū)別，所以最初錄制的語(yǔ)音信號(hào)不能滿足課設(shè)要求，后來(lái)使用3.1步驟中的錄音軟件解決了問(wèn)題。（2）由于對(duì)MATLAB使用能力有限，在語(yǔ)音信號(hào)加噪時(shí)出了錯(cuò)，頻譜圖中沒(méi)有出現(xiàn)單頻噪聲尖峰；然后通過(guò)網(wǎng)上搜索資料，并在同學(xué)的幫助下成功的顯示了正確的頻譜圖；在利用脈沖響應(yīng)不變法設(shè)計(jì)切比雪夫?yàn)V波器時(shí)得不到理想的濾波器，把通帶截止頻率與阻帶截止頻率之間的差值設(shè)置的太小或太大，通過(guò)適當(dāng)?shù)恼{(diào)整參數(shù)，繪制出來(lái)的圖像效果比較明顯，基本符合設(shè)計(jì)指標(biāo)。通過(guò)不斷調(diào)整As，最后取As=16dB,這樣通過(guò)MATLAB運(yùn)算出來(lái)的濾波器的阻帶波紋達(dá)到要求；在設(shè)計(jì)切比雪夫I型IIR濾波器過(guò)程中，程序運(yùn)行報(bào)錯(cuò)：Undefinedcommand/’afd_chbT，剛開(kāi)始不知道怎么回事兒，后來(lái)老師告訴我此錯(cuò)誤的原因在于沒(méi)有將自編調(diào)用函數(shù)放進(jìn)MATLAB的函數(shù)調(diào)用庫(kù)中，將調(diào)用函數(shù)’afd_chbT編寫(xiě)放入MATLAB的work文件夾后，問(wèn)題就解決了；5結(jié)束語(yǔ)課程設(shè)計(jì)的塵囂漸漸遠(yuǎn)去，懷著些許不舍，懷著對(duì)當(dāng)初課程設(shè)計(jì)開(kāi)始時(shí)候的豪情萬(wàn)丈的決心的留戀，懷著通過(guò)這次課程設(shè)計(jì)積累的信心與斗志，我開(kāi)始寫(xiě)下了這份結(jié)束語(yǔ)，為自己的足跡留下哪怕是微不足道但是對(duì)自己彌足珍貴的痕跡并期望與大家共勉。此次課程設(shè)計(jì)給我留下了太多記憶深刻的東西，學(xué)到了很多課本上學(xué)不到的知識(shí)。首先，我能夠運(yùn)用MATLAB，數(shù)字信號(hào)處理教程等學(xué)科所學(xué)的知識(shí)，并結(jié)合網(wǎng)絡(luò)，圖書(shū)館搜集到的各種各樣的材料，進(jìn)行分析，比較，取舍，對(duì)課本知識(shí)有了一個(gè)深入的了解，這對(duì)我的學(xué)業(yè)生涯是大有裨益的。其次，這次課程設(shè)計(jì)使我有了一次通過(guò)理論聯(lián)系實(shí)際，來(lái)解決實(shí)際問(wèn)題的經(jīng)歷，它培養(yǎng)了自己分析問(wèn)題，解決問(wèn)題的能力，以及上網(wǎng)檢索信息的能力；并且在和同學(xué)或者老師的交流中我體會(huì)到了團(tuán)隊(duì)的意義，以及探索未知學(xué)問(wèn)的樂(lè)趣。其實(shí)學(xué)到的知識(shí)其實(shí)是次要的，重要的是我們探索知識(shí)的過(guò)程，這個(gè)過(guò)程便是一個(gè)人自主學(xué)習(xí)能力的體現(xiàn)，它將影響著我今后的發(fā)展。最后，這次的課程設(shè)計(jì)也激發(fā)了我對(duì)通行工程本專業(yè)的興趣。曾經(jīng)我對(duì)我們通信工程專業(yè)是一片迷茫，但是通過(guò)這次課程設(shè)計(jì)，通過(guò)親自體驗(yàn)語(yǔ)音信號(hào)的加噪過(guò)濾過(guò)程，我似乎隱約了解了一些我們同形詞專業(yè)的未來(lái)的發(fā)展方向。這種感悟，這種興趣的來(lái)源只靠平時(shí)上課老師的說(shuō)教是肯定不可能領(lǐng)悟到的。參考文獻(xiàn)［1］程佩青.數(shù)字信號(hào)處理教程［M］.北京：清華大學(xué)出版社，2013⑵（美）維納K英格爾，（美）約翰G普羅克斯（著）.劉樹(shù)棠（譯）.數(shù)字信號(hào)處理（MATLAB版）.第2版.西安：西安交通大學(xué)出版社，2008［3］張威.MATLAB基礎(chǔ)與編程入門(mén)皿］.西安：西安電子科技大學(xué)出版社2008附錄1:語(yǔ)音信號(hào)濾波去噪設(shè)計(jì)源程序清單%程序名稱：finish.m%程序功能：語(yǔ)音信號(hào)去噪一一使用脈沖響應(yīng)不變法設(shè)計(jì)并聯(lián)型切比雪夫I型濾波器%程序作者：%最后修改日期：[x,fs,bits]=wavread('ldd.wav');%輸入?yún)?shù)為文件的全路徑和文件名，輸出的第一個(gè)參數(shù)是每個(gè)樣本的值，fs是生成該波形文件時(shí)的采樣率，bits是波形文件每樣本的編碼位數(shù)。sound(x,fs,bits);%按指定的采樣率和每樣本編碼位數(shù)回放N=length(x);%計(jì)算信號(hào)x的長(zhǎng)度f(wàn)n=2000;%單頻噪聲頻率，此參數(shù)可改t=0:1/fs:(N-1)/fs;%計(jì)算時(shí)間范圍，樣本數(shù)除以采樣頻率x=x';y=x+0.1*sin(fn*2*pi*t);sound(y,fs,bits);%應(yīng)該可以明顯聽(tīng)出有尖銳的單頻嘯叫聲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ì)算頻譜頻率范圍subplot(2,2,1);plot(t,x);title('原始信號(hào)');xlabel('時(shí)間/s');ylabel('幅度');axis([0,2.5,-1,1]);gridonsubplot(2,2,2);plot(f,X);title('語(yǔ)音信號(hào)幅度譜');xlabel('頻率，單位/Hz');ylabel('幅度譜');gridonsubplot(2,2,3);plot(t,y);title('加入一個(gè)單頻信號(hào)的語(yǔ)音信號(hào)');xlabel('時(shí)間，單位/s');ylabel('幅度');axis([0,2.5,-1.1,1.1]);gridonsubplot(2,2,4);plot(f,Y);title('加入一個(gè)單頻信號(hào)的語(yǔ)音信號(hào)幅度譜');xlabel('頻率，單位/Hz');ylabel('幅度譜');axis([0,4000,0,800]);gridonfs=8000;fn=3000;fp=fn-200;fc=fn-50;%定義通帶和阻帶截止頻率Rp=1;As=16;%定義通帶波紋和阻帶衰減wp=fp/fs*2*pi;ws=fc/fs*2*pi;%計(jì)算對(duì)應(yīng)的數(shù)字頻率T=1;%定義采樣間隔OmegaP=wp/T;OmegaS=ws/T;%截止頻率線性變換Ripple=10A(-Rp/20);Attn=10A(-As/20);%轉(zhuǎn)化為絕對(duì)形式[cs,ds]=afd_chb1(OmegaP,OmegaS,Rp,As);%計(jì)算濾波器階數(shù)和截止頻率[b,a]=imp_invr(cs,ds,T);%脈沖響應(yīng)不變法得到數(shù)字濾波器系統(tǒng)函數(shù)分子分母系數(shù)[db,mag,pha,grd,w]=freqz_m(b,a);%驗(yàn)證濾波器是否達(dá)到指定性能delta=[1,zeros(1,99)];ha=filter(b,a,delta);%計(jì)算脈沖響應(yīng)Subplot(221);plot(w/pi,db)%畫(huà)出濾波器的幅度響應(yīng)圖xlabel('w/pi');ylabel('dB');title('濾波器幅度響應(yīng)圖')Subplot(222);plot(w/pi,mag)%畫(huà)出濾波器幅度響應(yīng)圖xlabel('w/pi');ylabel('幅度mag');title('濾波器幅度響應(yīng)')Subplot(223);plot(w/pi,pha)%畫(huà)出濾波器的相位響應(yīng)圖xlabel('w/pi');ylabel('相位pha');title('濾波器相位響應(yīng)')Subplot(224);plot(ha)%畫(huà)出濾波器的脈沖響應(yīng)圖xlabel('n');ylabel('h(n)');title('濾波器脈