語音信號的數(shù)字濾波處理十一——橢圓函數(shù)hanning窗濾波器

1、 changsha university of science & technology數(shù)字信號處理課程設(shè)計 語音信號的數(shù)字濾波處理(十一) 題目： 橢圓函數(shù)（hanning窗）濾波器 學(xué)生姓名： 學(xué) 號： 班 級: 專 業(yè)： 指導(dǎo)教師： 實習(xí)起止時間: 2015年6月29日至2015年7月3日 語音信號的數(shù)字濾波處理(九) 題目 數(shù)字信號處理語音信號的數(shù)字濾波處理(十一)學(xué)生姓名： 學(xué) 號： 班 級： 所在院(系): 指導(dǎo)教師： 語音信號的數(shù)字濾波處理(十一) 摘要 本次設(shè)計的內(nèi)容為切比雪夫及hanning低通、高通、帶通濾波器，并利用matlab平臺進行設(shè)計。首先通過聲音處理語句得到聲音信

2、號的時域數(shù)據(jù)，利用fft變換可得到頻域數(shù)據(jù)，以此進行頻率分析。然后對原語音信號進行加噪處理，得到被污染信號。最后將被污染信號通過設(shè)計的濾波器，實現(xiàn)濾波功能，得到濾波后的語音信號。濾波器分別用切比雪夫ii型和hanning窗設(shè)計，間接法設(shè)計iir濾波器采用雙線性變換法，濾波器設(shè)計指標由頻譜分析得到。通過聲音播放語句進行語音播放，可觀察聲音的變化；通過圖形處理語句和fft得到時域圖和頻譜圖，可分析得到濾波器對頻率的濾波功能。關(guān)鍵詞：切比雪夫；聲音處理；hanning；matlab；fft目 錄1 緒論.21.1 課題背景.21.2 課題目的.22 課程設(shè)計預(yù)習(xí)與原理.32.1 課程設(shè)計預(yù)習(xí).32.

3、1.1 卷積運算的演示.32.2.2 采樣定理的演示.92.2 課程設(shè)計原理.112.2.1 頻譜分析原理.11 2.2.2 iir設(shè)計原理.12 2.2.3 fir設(shè)計原理.123 課程設(shè)計步驟和過程.153.1 iir設(shè)計步驟和過程.153.2 fir設(shè)計步驟和過程.154 設(shè)計程序的調(diào)試和運行結(jié)果.164.1切比雪夫低通濾波器程序的調(diào)試和運行結(jié)果.164.2切比雪夫高通濾波器程序的調(diào)試和運行結(jié)果.184.3切比雪夫帶通濾波器程序的調(diào)試和運行結(jié)果.204.4 hamming濾波器程序的調(diào)試和運行結(jié)果.235 總結(jié).25參考文獻.27附錄.28附錄 a.28附錄 b.33第 38 頁 共 5

4、6 頁1 緒論1.1 課題背景隨著軟硬件技術(shù)的發(fā)展，儀器的智能化與虛擬化已成為未來實驗室及研究機構(gòu)的發(fā)展方向。虛擬儀器技術(shù)的優(yōu)勢在于可由用戶定義自己的專用儀器系統(tǒng)，且功能靈活，很容易構(gòu)建，所以應(yīng)用面極為廣泛?；谟嬎銠C軟硬件平臺的虛擬儀器可代替?zhèn)鹘y(tǒng)的測量儀器，如示波器、邏輯分析儀、信號發(fā)生器、頻譜分析儀等。從發(fā)展史看，電子測量儀器經(jīng)歷了由模擬儀器、智能儀器到虛擬儀器，由于計算機性能的飛速發(fā)展，已把傳統(tǒng)儀器遠遠拋到后面，并給虛擬儀器生產(chǎn)廠家不斷帶來連鍋端的技術(shù)更新速率。目前已經(jīng)有許多較成熟的頻譜分析軟件，如spectralab、rsavu、dbfa等。聲卡是多媒體計算機最基本的配置硬件之一，價格

5、便宜，使用方便。matlab是一個數(shù)據(jù)分析和處理功能十分強大的工程實用軟件，他的數(shù)據(jù)采集工具箱為實現(xiàn)數(shù)據(jù)的輸入和輸出提供了十分方便的函數(shù)和命令。matlab是mathworks公司于1982年推出的一款功能強大、易于使用的高效數(shù)值計算和可視化軟件，它為進行算法開發(fā)、數(shù)據(jù)計算、信號分析與可視化提供了交互式應(yīng)用開發(fā)環(huán)境，主要包括基本數(shù)學(xué)計算、編程環(huán)境(m語言)、數(shù)據(jù)可視化、guide等。并附加了大量支持建模、分析、計算應(yīng)用的工具箱來擴展matlab基本環(huán)境，用于解決特定領(lǐng)域的工程問題，如今matlab已廣泛應(yīng)用于通信、信號處理、生物醫(yī)學(xué)、自動控制等領(lǐng)域。數(shù)字濾波器(digital filter,

6、df)在信號的過濾、檢測和參數(shù)估計等方面起著重要的作用。信號往往夾雜著噪聲及無用信號成分，必須將這些干擾成分濾除。數(shù)字濾波器對信號進行篩選，可通過特定頻段的信號。一般來說，噪聲信號往往是高頻信號，而經(jīng)典濾波器正是假定有用信號與噪聲信號具有不同的頻段，所以利用經(jīng)典濾波器可以去除噪聲。但如果有用信號和無用信號，或有用信號和噪聲的頻譜相互重疊，那么經(jīng)典濾波器則不能實現(xiàn)理想的濾波性能?，F(xiàn)代濾波器的作用是從含有噪聲的信號中估計出信號的某些特征或信號本身，一旦信號被估計出，那么估計出來的信號與原信號相比會有更高的信噪比。這類濾波器主要有維納濾波器、卡爾曼濾波器、線性預(yù)測濾波器以及自適應(yīng)濾波器等。1.2 課

7、題目的綜合運用數(shù)字信號處理的理論知識進行頻譜分析和濾波器設(shè)計，通過理論推導(dǎo)得出相應(yīng)結(jié)論，再利用 matlab 作為編程工具進行計算機實現(xiàn)，從而加深對所學(xué)知識的理解，建立概念。2.課程設(shè)計預(yù)習(xí)與原理2.1 課程設(shè)計預(yù)習(xí)2.1.1 卷積運算的演示(1) 線性卷積序列x1(n)=2 0 1 2 5 7 0 5 0 1 1 2，序列x2(n)=2 0 1 2 5 7 0 5。動態(tài)演示兩個序列進行線性卷積x1(n)x2(n)的翻轉(zhuǎn)、移位、乘積、求和的過程。其中翻轉(zhuǎn)采用fliplr1，程序如下：n=-7:18; m=17;yn=zeros(1,19);figure(1)stem(yn);xlabel(n)

8、ylabel(y(n)xn1=2 0 1 2 5 7 0 5 0 1 1 2;xm1=zeros(1,7) xn1 zeros(1,7);figure(2)stem(n,xm1)xlabel(m)ylabel(x1(m)xn2=2 0 1 2 5 7 0 5;xm2=fliplr(xn2) zeros(1,18);figure(3)stem(n,xm2)xlabel(m)ylabel(x2(n-m)title(n=0)yn(1)=sum(xm1.*xm2);figure(4)stem(yn)xlabel(n)ylabel(y(n)title(n=n)for n=1:17xm3=zeros(1,

9、n) fliplr(xn2) zeros(1,m);figure(5)stem(n,xm3)xlabel(m)ylabel(x2(n-m)title(n=n)m=m-1;yn(n+1)=sum(xm1.*xm3);figure(6)stem(yn)xlabel(n)ylabel(y(n)title(n=n)endxm3=zeros(1,18) fliplr(xn2)figure()stem(xm3);xlabel(m)ylabel(x2(n-m)title(n=n)yn(19)=sum(xm1.*xm3);figure()stem(yn)xlabel(n)ylabel(y(n)title(n=

10、n)線性卷積動態(tài)演示過程見：附錄a。(2) 循環(huán)卷積序列x1(n)=2 0 0 7 5 7 1 7 0 1 1 9，序列x2(n)=2 0 0 7 4 7 0 3，n=12。動態(tài)演示兩個序列進行圓周卷積2, 3x1(n)x2(n)的翻轉(zhuǎn)、移位、乘積、求和的過程。程序如下：n=0:11; yn=zeros(1,12);figure(1)stem(yn);xlabel(n)ylabel(y(n)xn1=2 0 0 7 5 7 1 7 0 1 1 9;figure(2)stem(n,xn1)xlabel(m)ylabel(x1(m)xn2=2 0 0 7 4 7 0 3;xm2=xn2 zeros(

11、1,length(xn1)-length(xn2);figure(3)stem(n,xm2)xlabel(m)ylabel(x2(m)title(n=0);yn(1)=sum(xn1.*xm2);figure(4)stem(yn)xlabel(n)ylabel(y(n)title(n=n)for n=1:11xm1=fliplr(xn1(1:n) fliplr(xn1(n+1:12);figure()stem(n,xm1)xlabel(m)ylabel(x1(n-m)title(n=n)yn(n)=sum(xm1.*xm2);figure()stem(n,yn)xlabel(n)ylabel(

12、y(n)title(n=n)endfigure()xm1=fliplr(xn1);stem(n,xm1)xlabel(m)ylabel(x1(n-m)title(n=n)yn(12)=sum(xm1.*xm2);figure()stem(n,yn)xlabel(n)ylabel(y(n)title(n=n)當n12時，圓周卷積發(fā)生混疊；當n20時，圓周卷積等于線性卷積。圓周卷積動態(tài)演示過程見：附錄a。(3) 聲音文件線性卷積序列x1(n)=2 0 1 2 5 7 0 5,讀取一段聲音數(shù)據(jù)（windows xp 菜單命令），當循環(huán)卷積長度大于或等于兩序列長度之和時，循環(huán)卷積等于線性卷積。因為直接

13、用fft進行1024點卷積大于兩序列長度，所以可用線性卷積替代圓周卷積，其程序如下：n=-7:686; m=685;yn=zeros(1,687);xn1,fs,nbits=wavread(g:windows xp 菜單命令);xn1(:,2)=xm1=zeros(1,7) xn1 zeros(1,7);xn2=2 0 1 2 5 7 0 5;xm2=fliplr(xn2) zeros(1,686);xm1=xm1(:,1:694)yn(1)=sum(xm1.*xm2);for n=1:685xm3=zeros(1,n) fliplr(xn2) zeros(1,m);m=m-1;yn(n+1)

14、=sum(xm1.*xm3)endxm3=zeros(1,686) fliplr(xn2);yn(687)=sum(xm1.*xm3);figure(1)stem(yn)xlabel(m)ylabel(y(n)title(n=n)線性卷積結(jié)果如圖2.1所示。 圖2.1 線性卷積結(jié)果(4) fft 由循環(huán)卷積定理可知：對于時域序列循環(huán)卷積，可先進行fft變換，然后頻率相乘，最后對結(jié)果進行ifft變換，即可得到時域循環(huán)卷積結(jié)果。其程序如下：程序：y,fs,nbits=wavread(g:windows xp 菜單命令);y=fft(y,1024);xn2=2 0 1 2 5 7 0 5;x=fft

15、(xn2,1024);x1=rot90(x,3);z=x1.*y;z=ifft(z,1024)figure(2)stem(z)axis(0,700,-1,0.6);fft求得結(jié)果如圖2.2所示。 圖2.2 1024點fft求得結(jié)果由圖2.1和圖2.2可得：直接用卷積求得結(jié)果與間接用fft求得結(jié)果相同。同時可觀察到間接用fft求得速度要比直接用卷積求得速率要快，效率要高。2.1.2采樣定理的演示 編寫程序演示采樣定理（時域采樣、頻譜周期延拓），同時演示采樣頻率小于2fc時，產(chǎn)生的混疊效應(yīng)： 對下面連續(xù)信號進行采樣：，a為幅度因子，a為衰減因子，為模擬角頻率，其中n為學(xué)號(例如，王墨同學(xué)n=23)

16、程序：fs=25;n=12;t=0:1/fs:1;t1=0:0.0005:1;a=n;a=sqrt(2)*n*pi;w=a;x=a*exp(-a*t).*sin(w*t);x1=a*exp(-a*t1).*sin(w*t1);subplot(5,1,1)plot(x1);xn=ones(1,length(t);subplot(512)stem(xn)yn=xn.*x;subplot(5,1,3)stem(yn);y=fft(yn,fs);f=fs*(0:fs/2)/fs;subplot(5,1,4)plot(f,y(1:fs/2+1);yn=ifft(y);subplot(5,1,5)plot

17、(yn)其結(jié)果如果2.3.1所示；圖2.3.1 采樣定理演示1 要求輸入采樣頻率fs（根據(jù)程序處理需要指定范圍）后，在時域演示信號波形、采樣脈沖及采樣后信號；在頻域演示不同采樣頻率下對應(yīng)信號的頻譜。程序：fs=60n=12t=0:1/fs:1;t1=0:0.0005:1;a=n;a=sqrt(2)*n*pi;w=a;x=a*exp(-a*t).*sin(w*t);x1=a*exp(-a*t1).*sin(w*t1);subplot(5,1,1)plot(x1);xn=ones(1,length(t);subplot(512)stem(xn)yn=xn.*x;stem(yn);y=fft(yn,

18、fs);f=fs*(0:fs/2)/fs;subplot(5,1,4)plot(f,y(1:fs/2+1);yn=ifft(y);subplot(5,1,5)plot(yn)其結(jié)果如果2.3.2所示；. 圖2.3.2 采樣定理演示22.2 課程設(shè)計原理2.2.1 頻譜分析原理時域分析只能反映信號的幅值隨時間的變化情況，除單頻率分量的簡單波形外，很難明確提示信號的頻率組成和各頻率分量大小，而頻譜分析能很好的解決此問題。(1) dft與fft對于給定的時域信號y，可以通過fourier變換得到頻域信息y。y可按下式計算 （2-1）式中，n為樣本容量，為采樣間隔。采樣信號的頻譜是一個連續(xù)的頻譜，不可

19、能計算出所有的點的值，故采用離散fourier變換(dft)，即 （2-2）式中，。但上式的計算效率很低，因為有大量的指數(shù)(等價于三角函數(shù))運算，故實際中多采用快速fourier變換(fft)。其原理即是將重復(fù)的三角函數(shù)算計的中間結(jié)果保存起來，以減少重復(fù)三角函數(shù)計算帶來的時間浪費。由于三角函數(shù)計算的重復(fù)量相當大，故fft能極大地提高運算效率。(2) 頻譜圖為了直觀地表示信號的頻率特性，工程上常常將fourier變換的結(jié)果用圖形的方式表示，即頻譜圖。以頻率f為橫坐標，為縱坐標，可以得到幅值譜；以頻率f為橫坐標，為縱坐標，可以得到相位譜；以頻率f為橫坐標，為縱坐標，可以得到實頻譜；以頻率f為橫坐標

20、，為縱坐標，可以得到虛頻譜。根據(jù)采樣定理，只有頻率不超過的信號才能被正確采集，即fourier變換的結(jié)果中頻率大于的部分是不正確的部分，故不在頻譜圖中顯示。即橫坐標。2.2.2 iir設(shè)計原理(1) 切比雪夫濾波器原理切比雪夫濾波器的幅頻特性具有等波紋特性。它有兩種形式：振幅特性在通帶內(nèi)是等波紋的，在阻帶內(nèi)是單調(diào)下降的切比雪夫i型；振幅特性在阻帶內(nèi)是等波紋的，在通帶內(nèi)是單調(diào)下降的切比雪夫ii型。(2) 雙線性變換法工作原理 雙線性變換中數(shù)字域頻率和模擬頻率之間的非線性關(guān)系限制了它的應(yīng)用范圍，只有當非線性失真是允許的或能被忽略時，才能采用雙線性變換法，通常低通、高通、帶通和帶阻等濾波器等具有分段

21、恒定的頻率特性，可以采用預(yù)畸變的方法來補償頻率畸變，因此可以采用雙線性變換設(shè)計方法。(3) 脈沖響應(yīng)不變法工作原理 沖激響應(yīng)不變法遵循的準則是使數(shù)字濾波器的單位取樣響應(yīng)與參照的模擬濾波器的脈沖響應(yīng)的取樣值完全一樣，即h(n)=ha(nt),其中t為取樣周期。實際是由模擬濾波器轉(zhuǎn)換成為數(shù)字濾波器，就是要建立模擬系統(tǒng)函數(shù)ha（s）與數(shù)字系統(tǒng)函數(shù)h（z）之間的關(guān)系。脈沖響應(yīng)不變法是從s平面映射到z平面，這種映射不是簡單的代數(shù)映射，而是s平面的每一條寬為2/t的橫帶重復(fù)地映射到整個z平面。 2.2.3 fir設(shè)計原理由于iir數(shù)字濾波器能夠保留一些模擬濾波器的優(yōu)良特性，因此應(yīng)用很廣。但是這些特性是以犧

22、牲線性相位頻率特性為代價的，即用butterworth、切比雪夫和橢圓法設(shè)計的數(shù)字濾波器逼近理想的濾波器的幅度頻率特性，得到的濾波器往往是非線性的。在許多電子系統(tǒng)中，對幅度頻率特性和線性相位特性都有較高的要求，所以iir濾波器在這些系統(tǒng)中往往難以勝任。有限長單位沖激響應(yīng)(fir)數(shù)字濾波器具有以下優(yōu)良的特點：(1) 可在設(shè)計任意幅度頻率特性濾波器的同時，保證精確、嚴格的線性相位特性。(2) fir數(shù)字濾波器的單位沖激響應(yīng)h(n)是有限長的，可以用一個固定的系統(tǒng)來實現(xiàn)，因而fir數(shù)字濾波器可以做成因果穩(wěn)定系統(tǒng)。 (3) 允許設(shè)計多通帶(多阻帶)系統(tǒng)。 窗函數(shù)法就是設(shè)計fir數(shù)字濾波器的最簡單的方

23、法。它在設(shè)計fir數(shù)字濾波器中有很重要的作用，正確地選擇窗函數(shù)可以提高設(shè)計數(shù)字濾波器的性能，或者在滿足設(shè)計要求的情況下，減小fir數(shù)字濾波器的階次。常用的窗函數(shù)4有以下幾種：矩形窗(rectangular window)、三角窗(triangular window)、漢寧窗(hanning window)、海明窗(hamming window)、布拉克曼窗(blackman window)、切比雪夫窗(chebyshev window)、巴特里特窗(bartlett window)及凱塞窗(kaiser window)。在matlab中，實現(xiàn)矩形窗的函數(shù)為boxcar和rectwin，其調(diào)用格

24、式如下： w=boxcar(n)w=rectwin(n)其中n是窗函數(shù)的長度，返回值w是一個n階的向量，它的元素由窗函數(shù)的值組成。實際上，w=boxcar(n)等價于w=ones(n,1)。在matlab中，實現(xiàn)三角窗的函數(shù)為triang，調(diào)用格式為： w=triang(n) 在matlab中，實現(xiàn)漢寧窗的函數(shù)為hann，調(diào)用格式如下： w=hann(n)w=hann(n,sflag)hann函數(shù)中的參數(shù)sflag為采樣方式，其值可取symmetric(默認值)或periodic。當sflagsymmetric時，為對稱采樣；當sflagperiodic時，為周期采樣，此時hann函數(shù)計算n+

25、1個點的窗，但是僅返回前n個點。在matlab中，實現(xiàn)海明窗的函數(shù)為hamming，調(diào)用格式分別如下： w=hamming (n)w=hamming (n,sflag)其中sflag的用法同上。在matlab中，實現(xiàn)布拉克曼窗的函數(shù)為blackman，調(diào)用格式如下： w=blackman (n)w=blackman (n,sflag)在matlab中，實現(xiàn)切比雪夫窗的函數(shù)為chebwin，調(diào)用格式為： w=chebwin (n,r)其中r 表示切比雪夫窗函數(shù)的傅里葉變換旁瓣幅度比主瓣低rdb(其默認值為100db)，且旁瓣是等紋波的。在matlab中，實現(xiàn)巴特里特窗的函數(shù)為bartlett，調(diào)

26、用格式為： w=bartlett (n)在matlab中，實現(xiàn)凱塞窗的函數(shù)為kaiser，調(diào)用格式為： w=kaiser (n,beta)其中beta為窗函數(shù)的參數(shù)。各種窗函數(shù)的性能比較可見表2-1，在設(shè)計fir濾波器的過程中可以根據(jù)要求選擇合適的窗函數(shù)：表2-1 各種窗函數(shù)的性能比較窗 函 數(shù)第一旁瓣相對于主瓣衰減/db主 瓣 寬阻帶最小衰減/db矩形窗134/n21三角窗258/n25漢寧窗318/n44海明窗418/n53布拉克曼窗5712/n74凱塞窗可調(diào)可調(diào)可調(diào)切比雪夫窗可調(diào)可調(diào)可調(diào)3 課程設(shè)計步驟和過程3.1 iir設(shè)計步驟和過程(1) 按一定的規(guī)則將給出的數(shù)字濾波器的技術(shù)指標轉(zhuǎn)換

27、為模擬低通濾波器的技術(shù)指標。(2) 根據(jù)轉(zhuǎn)換后的技術(shù)指標設(shè)計模擬低通濾波器g(s)。(3) 再按一定的規(guī)則將g(s)轉(zhuǎn)換成h(s)。若所設(shè)計的數(shù)字濾波器是低通的，那么上述設(shè)計工作可以結(jié)束，若所設(shè)計是高通、帶通或帶阻濾波器，那么還需進行以下步驟。(4) 將高通、帶通或帶阻數(shù)字濾波器的技術(shù)指標先轉(zhuǎn)換為低通模擬濾波器的技術(shù)指標，然后按照上述步驟(2)設(shè)計出低通g(s)，再將g(s)轉(zhuǎn)換為所需的h(z)。由于的頻率映射關(guān)系是根據(jù)推導(dǎo)的，所以使j軸每隔2/ts便映射到單位圓上一周，利用沖激響應(yīng)不變法設(shè)計數(shù)字濾波器時可能會導(dǎo)致上述的頻域混疊現(xiàn)象。為了克服這一問題，需要找到由s平面到z平面的另外的映射關(guān)系，

28、這種關(guān)系應(yīng)保證：(1) s平面的整個j軸僅映射為z平面單位圓上的一周；(2) 若g(s)是穩(wěn)定的，由g(s)映射得到的h(z)也應(yīng)該是穩(wěn)定的；(3) 這種映射是可逆的，既能由g(s)得到h(z)，也能由h(z)得到g(s)；(4) 如果g(j0)=1，那么。雙線性z變換滿足以上4個條件的映射關(guān)系，其變換公式為 (3-1)雙線性z變換的基本思路是：首先將整個s平面壓縮到一條從/ts到/ts的帶寬為2/ts的橫帶里，然后通過標準的變換關(guān)系將橫帶變換成整個z平面上去，這樣就得到s平面與z平面間的一一對應(yīng)的單值關(guān)系。3.2 fir課程設(shè)計步驟和過程(1)根據(jù)對阻帶衰減及過渡帶的指標要求，選擇窗函數(shù)的類

29、型，并估計窗口長度n；(2)構(gòu)造希望逼近的頻率響應(yīng)函數(shù)hd(ejw)，即hd(ejw)=hdg(w)e-jw（n-1）/2；(3)計算hd(n)；(4)加窗得到設(shè)計結(jié)果：h(n)=hd(n)w(n)。數(shù)字高通濾波器的設(shè)計分兩種情況，即無相移和相移為/2兩種情況。當無相移時，n為奇數(shù)時，所設(shè)計的fir數(shù)字高通濾波器為型濾波器；當n為偶數(shù)時，為型濾波器。當相移為/2時，n為奇數(shù)時，所設(shè)計的fir數(shù)字高通濾波器為型濾波器；當n為偶數(shù)時，為型濾波器。數(shù)字帶通濾波器的設(shè)計也分兩種情況討論，即無相移和相移為/2兩種情況。當無相移時，n為奇數(shù)時，所設(shè)計的fir數(shù)字帶通濾波器為型濾波器；當n為偶數(shù)時，為型濾波

30、器。當相移為/2時，n為奇數(shù)時，所設(shè)計的fir數(shù)字帶通濾波器為型濾波器；當n為偶數(shù)時，為型濾波器。4 設(shè)計程序的調(diào)試和運行結(jié)果4.1 切比雪夫低通程序的調(diào)試和運行結(jié)果 如圖4.1所示為原語音信號的時域圖和頻譜圖。在matlab中通過wavplay(filename)讀取語音信號數(shù)據(jù)。其程序如下：y,fs,nbits=wavread(g:windows background.wav); 圖4.1 原語音信號時域圖和頻譜圖如圖4.2所示為被污染語音信號的時域圖和頻譜圖。其加噪信號程序如下：t=0:0.000045:0.90855;x1=0.06*sin(2*pi*5500*t)+0.07*sin(

31、2*pi*6000*t); 圖4.2 被污染語音信號時域圖和頻譜圖如圖4.3所示為低通濾波器的幅頻圖。設(shè)計指標：wpz=4500hz; wsz=4700hz; rp=2; rs=100;其設(shè)計程序如下：wpz=2*4500/fs;wsz=2*4700/fs;wp=2*fs*tan(wpz*pi/2);ws=2*fs*tan(wsz*pi/2);rp=2;rs=100;n,wc=cheb2ord(wp,ws,rp,rs,s);b,a=cheby2(n,rs,wc,s);bz,az=bilinear(b,a,fs);hk=freqz(bz,az,16384,fs); 圖4.3 低通濾波器幅頻圖如圖

32、4.4所示為濾波后語音信號的時域圖和頻譜圖。其濾波程序如下：yk=z.*hk; 圖4.4 濾波后語音信號的時域圖和頻譜圖4.2 切比雪夫高通程序的調(diào)試和運行結(jié)果 如圖4.5所示為原語音信號的時域圖和頻譜圖。在matlab中通過wavplay(filename)讀取語音信號數(shù)據(jù)。其程序如下：y,fs,nbits=wavread(g:windows background.wav); 圖4.5 原語音信號的時域圖和頻譜圖如圖4.6所示為被污染語音信號的時域圖和頻譜圖。其加噪信號程序如下：t=0:0.000045:0.90855;x1=0.1*sin(2*pi*100*t)+0.2*sin(2*pi*

33、220*t);x=rot90(x1,3),rot90(x1,3);z=x+y;圖4.6 被污染語音信號的時域圖和頻譜圖如圖4.7所示為高通濾波器的幅頻圖。設(shè)計指標：wpz=700hz; wsz=500hz; rp=10; rs=100;其設(shè)計程序如下：wpz=2*700/fs;wsz=2*500/fs;wp=2*fs*tan(wpz*pi/2);ws=2*fs*tan(wsz*pi/2);rp=10;rs=100;n,wc=cheb2ord(wp,ws,rp,rs,s);b,a=cheby2(n,rs,wc,high,s);bz,az=bilinear(b,a,fs);hk=freqz(bz,

34、az,16384,fs);圖4.7 高通濾波器的幅頻圖如圖4.8所示為濾波后語音信號的時域圖和頻譜圖。其濾波程序如下：yk=z.*hk;圖4.8 濾波后語音信號的時域圖和頻譜圖4.3 切比雪夫帶通程序的調(diào)試和運行結(jié)果如圖4.9所示為原語音信號的時域圖和頻譜圖。在matlab中通過wavplay(filename)讀取語音信號數(shù)據(jù)。其程序如下：y,fs,nbits=wavread(g:windows background.wav);圖4.9 原語音信號的時域圖和頻譜圖如圖4.10所示為被污染語音信號的時域圖和頻譜圖。其加噪信號程序如下：t=0:0.000045:0.90855;x1=0.1*si

35、n(2*pi*100*t)+0.2*sin(2*pi*220*t)+ 0.07*sin(2*pi*6000*t);x=rot90(x1,3),rot90(x1,3);z=x+y; 圖4.10 被污染語音信號的時域圖和頻譜圖 如圖4.11所示為帶通濾波器的幅頻圖。設(shè)計指標：wpz=4500hz; wsz=4700hz; wpzl=700hz; wszl=500hz; rp=20; rs=50;其設(shè)計程序如下：wpz=2*700/fs,2*4500/fs;wsz=2*500/fs,2*4700/fs;wp=2*fs*tan(wpz*pi/2);ws=2*fs*tan(wsz*pi/2);rp=20

36、;rs=50;n,wc=cheb2ord(wp,ws,rp,rs,s);b,a=cheby2(n,rs,wc,s);bz,az=bilinear(b,a,fs);hk=freqz(bz,az,16384,fs);圖4.11 帶通濾波器的幅頻圖如圖4.12所示為濾波后語音信號的時域圖和頻譜圖。其濾波程序如下：yk=z.*hk;圖4.12 濾波后語音信號的時域圖和頻譜圖4.4 hanning程序的調(diào)試和運行結(jié)果如圖4.13所示為原語音信號的時域圖和頻譜圖。在matlab中通過wavplay(filename)讀取語音信號數(shù)據(jù)。其程序如下：y,fs,nbits=wavread(g:windows b

37、ackground.wav);圖4.4.1 原語音信號的時域圖和頻譜圖圖4.4.2 帶通濾波器的幅頻圖和相頻圖圖4.4.3 低通濾波器的幅頻圖和相頻圖圖4.4.4高通濾波器的幅頻圖和相頻圖5 總結(jié)在這次課程設(shè)計中，我遇到了許多問題，這些問題是我們初學(xué)者經(jīng)常會遇到的問題。通過老師的指導(dǎo)，和自己的思考，總結(jié)如下：問題1：subscript indices must either be real positive integers or logicals.解決辦法：在訪問矩陣（包括向量、二維矩陣、多維數(shù)組，下同）的過程中，下標索引要么從 0 開始，要么出現(xiàn)了負數(shù)。注：matlab 的語法規(guī)定矩陣的索

38、引從 1 開始，這與 c 等編程語言的習(xí)慣不一樣。調(diào)試一下程序，把下標為 0 或者負數(shù)的地方修正。問題2：undefined function or variable a.解決辦法：1. 如果 a 是變量，則檢查變量是否被賦值。2. 如果 a 是函數(shù)，則這是因為matlab在所有已添加的路徑中都無法找到該函數(shù)對應(yīng)的m文件而導(dǎo)致的，對此，把該m文件移動到當前路徑下再運行便可7。問題3：inner matrix dimensions must agree，由于各種原因使得矩陣長度不匹配，即使長度一樣，但因為一個是行矩陣，一個是列矩陣，程序報錯。解決方法：長度不一致時，將短的矩陣增零，使長度一致；矩

39、陣類型不同時，使用rot90(xn,3)：矩陣逆時旋轉(zhuǎn)270度，使行轉(zhuǎn)換成列或列轉(zhuǎn)換成行，或者直接使用xn直接轉(zhuǎn)換。問題4：經(jīng)過很多次檢查，程序沒錯，有時候輸入課本上的程序也會報錯。解決方法：檢查打開文件目錄，一般打開文件需要在matble下的work內(nèi)。matlab中常用語句：(1) 學(xué)會用help和doc函數(shù)(2) 輸入輸出文件：save/load，在屏幕上顯示文件：type(3) 解線性方程組ax=b：x=ab(4) 作圖時兩張曲線合并：hold on或者subplot作子圖(5) 程序計算時間：tic，toc或者clock(6) 變量顯示方式更改：format long/short/b

40、ank(7) 數(shù)組元素求和：sum(8) 求數(shù)組長度：length，求矩陣維數(shù)：size或者ndims，矩陣元素個數(shù)：numel(9) 函數(shù)作圖：餅圖：pie/pie3，誤差圖：errorbar，散點圖：scatter/scatter3，直方圖：hist，函數(shù)圖：fplot，動畫：movie(10) 矩陣分析：左右翻轉(zhuǎn)：fliplr，上下翻轉(zhuǎn)：flipud，轉(zhuǎn)置：transpose，矩陣求逆：inv，矩陣范數(shù)：norm，條件數(shù)：cond，初等變換：rref，特征值：eig/eigs(11) 特殊矩陣：元素全為1的矩陣：ones，元素全為0的矩陣：zeros，單位陣：eye，魔方陣：magic，

41、線性變化數(shù)組：linspace，聚合矩陣：cat/horzcat/vertcat8(12) 隨機數(shù)：創(chuàng)建一個元素服從均勻分布的隨機數(shù)數(shù)組：rand，創(chuàng)建一個元素服從正態(tài)分布的隨機數(shù)數(shù)組：randn，二項分布：binornd 指數(shù)分布：exprnd， f分布：frnd，幾何分布：geornd，超幾何分布：hygernd，泊松分布：poissrnd，正態(tài)分布：normrnd，離散均勻分布：unidrnd，連續(xù)均勻分布：unifrnd(13) 清屏：clc 清理內(nèi)存：clear(14) 語句太長的話可以再句末加.換行 通過近兩周的課程設(shè)計，我學(xué)到了許多知識，使我實現(xiàn)了從茫然到基本的認識到比較熟練的運

42、用matlab的轉(zhuǎn)變。剛開始自己編的程序運行就錯誤很多，當時都沒信心能把這個課程設(shè)計完成，但什么都要自己面對的。后來在老師的指導(dǎo)和自己的思考下，我把程序統(tǒng)統(tǒng)都重新寫了一遍，并反復(fù)修改，最后通過編譯。在后面的時間里，我順利完成了濾波器設(shè)計的功能。參考文獻1 徐東艷，孟曉剛matlab函數(shù)庫查詢辭典m北京：中國鐵道出版社，20062 劉衛(wèi)國matlab程序設(shè)計應(yīng)用m北京：高等教育出版社，20083 趙紅怡數(shù)字信號處理及matlab實現(xiàn)m北京：化學(xué)工業(yè)出版社，20024 高西全，丁玉美數(shù)字信號處理（第三版）m陜西：西安電子科技大學(xué)出版社，20105 何強，何英matlab擴展編程m北京：清華大學(xué)出版

43、社，20026 車晴電子系統(tǒng)仿真與matlabm北京：北京廣播學(xué)院出版社，20027 孫洪，余翔宇數(shù)字信號處理實驗指導(dǎo)書m北京：電子工業(yè)出版社，20078 vinay k.ingle，john g.proakisdigital signal processing using matlabm北京：科學(xué)出版社，2006附錄a 線性卷積和圓周卷積動態(tài)演示過程（1）線性卷積動態(tài)演示 圖a25 x2(3-m) 圖a26 x2(4-m) 圖a27 x2(5-m) 圖a28 x2(6-m) 圖a29 x2(7-m) 圖a30 x2(8-m) 附錄b 濾波器完整程序(1) 切比雪夫低通濾波器完整程序：y,fs

44、,nbits=wavread(g:windows background.wav);figure(1)subplot(2,1,1);plot(abs(y);xlabel(n);ylabel(y);title(原語音信號時域波形)y=fft(y,32768);f=22050*(0:16384)/32768;subplot(2,1,2);plot(f,abs(y(1:16385);xlabel(f);ylabel(y);title(原語音信號頻譜圖)t=0:0.000045:0.90855;x1=0.06*sin(2*pi*5500*t)+0.07*sin(2*pi*6000*t);x=rot90(x1,3),rot90(x1,3);y=y(1:20191,:)z=x+y;f