1、 語音信號的數字濾波處理(九) 1 緒論1.1 課題背景數字濾波是數字信號處理的重要內容，數字濾波器可分為IIR和FIR兩大類。對于IIR數字濾波器的設計，需要借助模擬原型濾波器，再將模擬濾波器轉化為數字濾波器，文中采用的設計方法是脈沖響應不變法、雙向性變換法和完全函數設計法；對于FIR數字濾波器的設計，可以根據所給定的頻率特性直接設計，文中采用的設計方法是窗函數法。根據IIR濾波器和FIR濾波器的特點，在MATLAB壞境下分別用雙線性變換法設計IIR和用窗函數設計FIR數字濾波器，并對采集的語音信號進行分析，最后給出了IIR和FIR對語音濾波的效果。1.2 課題要求1. 掌握數字信號處理的基

2、本概念，基本理論和基本方法。2. 熟悉離散信號和系統的時域特性。3. 掌握序列快速傅里葉變換方法。4. 學會MATLAB的使用，掌握MATLAB的程序設計方法。5. 掌握利用MATLAB對語音信號進行頻譜分析。6. 掌握濾波器的網絡結構。 2 課程設計預習與原理 2.1 卷積運算的演示 2.1.1 線性卷積 序列x1(n)=2 0 1 2 5 7 0 5 0 2 0 3，序列x2(n)=2 0 0 1 0 1 1 0。動態演示兩個序列進行線性卷積x1(n)x2(n)的翻轉、移位、乘積、求和的過程。其中翻轉采用fliplr，程序如下：n=-7:18;M=17;yn=zeros(1,19);fig

3、ure(1) stem(yn);xlabel(n)ylabel(y(n)xn1=2 0 1 2 5 7 0 5 0 2 0 3;xm1=zeros(1,7) xn1 zeros(1,7);%為26個數字的矩陣figure(2)stem(n,xm1)xlabel(m)ylabel(x1(m)xn2=2 0 0 1 0 1 1 0; xm2=fliplr(xn2) zeros(1,18); %移位，補零為26個數字的矩陣figure(3)stem(n,xm2)xlabel(m)ylabel(x2(N-m)title(n=0)yn(1)=sum(xm1.*xm2);%對xm1與xm2進行對應原素乘方

4、之后進行數組轉置，求和;即為求卷積figure(4)stem(yn)xlabel(n)ylabel(y(n)title(n=N)for N=1:17 xm3=zeros(1,N) fliplr(xn2) zeros(1,M); figure();stem(n,xm3)xlabel(m)ylabel(x2(N-m)title(n=N)M=M-1;yn(N+1)=sum(xm1.*xm3);figure()stem(yn)xlabel(n)ylabel(y(n)title(n=N)end xm3=zeros(1,18) fliplr(xn2)figure()stem(xm3)xlabel(m)yl

5、abel(x2(N-m)title(n=N);yn(19)=sum(xm1.*xm3);figure()stem(yn)xlabel(n)ylabel(y(n) 2.1.2循環卷積 序列x1(n)=2 0 1 2 5 7 0 5 0 2 0 3，序列x2(n)=2 0 0 1 0 1 1 0，N=12。動態演示兩個序列進行圓周卷積x1(n)x2(n)的翻轉、移位、乘積、求和的過程。程序如下：n=0:11;yn=zeros(1,12);figure(1)stem(yn)xlabel(n)ylabel(y(n) %圖1,12個0 xn1=2 0 1 2 5 7 0 5 0 2 0 3; figur

6、e(2)stem(n,xn1)xlabel(m)ylabel(x1(m) xn2=2 0 0 1 0 1 1 0; xm2=xn2 zeros(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)yla

7、bel(x1(N-m)title(n=N)yn(N)=sum(xm1.*xm2);figure()stem(n,yn)xlabel(n)ylabel(y(n)title(n=N)end figure()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) 2.1.3聲音文件線性卷積 序列x1(n)=2 0 0 1 0 1 1 0,讀取一段聲音數據，當循環卷積長度大于或等于兩序列長度之和時，

8、循環卷積等于線性卷積。因為直接用FFT進行1024點卷積大于兩序列長度，所以可用線性卷積替代圓周卷積，其程序如下：n=0:686;M=2053;yn=zeros(1,687);xn1,fs,nbits=wavread(C:UsersacerDesktop鄧紫棋-你把我灌醉.wav);xn1=xn1(:);xn1=xn1;xn1=xn1(10000:12047);%xn1=xn1;xm1=zeros(1,7) xn1 zeros(1,7);xn2=2 0 0 1 0 1 1 0;xm2=fliplr(xn2) zeros(1,2054);yn(1)=sum(xm1.*xm2);for N=1:6

9、85xm3=zeros(1,N) fliplr(xn2) zeros(1,M);M=M-1;yn(N+1)=sum(xm1.*xm3)endxm3=zeros(1,2054) fliplr(xn2);yn(687)=sum(xm1.*xm3);figure(1)stem(yn)xlabel(m)ylabel(y(n)title(n=N)線性卷積結果如圖所示 圖1 由循環卷積定理可知：對于時域序列循環卷積，可先進行FFT變換，然后頻率相乘，最后對結果進行IFFT變換，即可得到時域循環卷積結果。其程序如下：y,fs,nbits=wavread(C:UsersacerDesktop鄧紫棋-你把我灌醉

10、.wav);y=y(10000:12047);Y=fft(y,1024);xn2=1 2.43 6.17 12.93 22.17 32.25 40.88 45.87 45.87 40.88 32.25 22.17 12.93 6.17 2.43 1;X=fft(xn2,1024);X1=rot90(X,3);Z=(X1.*Y);z=ifft(Z,1024)*8figure(2)stem(z)axis(0,700,-1,0.6);結果如圖所示： 圖2 2.1.4采樣定理的演示 序列長度為28，先計算x(n)的512點FFT，得到其頻譜函數X512(k)，再對X512(k)隔點抽取，得到X32(k

11、)和X16(k)。再分別對其進行IFFT變換。觀察其時域圖。其程序如下：fs=27;n=3;t=0:1/fs:1;xa=n*exp(-n*sqrt(2)*pi)*sin(sqrt(2)*n*t);subplot(3,2,1)stem(xa)title(a) 原28點時域信號xa)Xk=fft(xa,512);subplot(3,2,2)stem(Xk)title(b) 512點FFTx(t)X32k=Xk(1:16:512);subplot(3,2,4)stem(X32k)title(d) 頻率采樣頻率為32時的頻譜圖)x32a=ifft(X32k);subplot(3,2,3)stem(x3

12、2a);title(c) 32點IFFTX32(k)得到的x32a(t)X16k=X32k(1:2:32);subplot(3,2,6)stem(X16k)title(f) 頻率采樣頻率為16時的頻譜圖)x16a=ifft(X16k);subplot(3,2,5)stem(x16a);title(e) 16點IFFTX16(k)得到的x16a(t)結果如圖所示： 圖3由圖可知：當進行32點IFFT時，無時域混疊失真；當進行16點IFFT時，有時域混疊失真。2.2 課程設計原理2.2.1 IIR設計原理(1) 切比雪夫濾波器原理切比雪夫濾波器的幅頻特性具有等波紋特性。它有兩種形式：振幅特性在通帶

13、內是等波紋的，在阻帶內是單調下降的切比雪夫I型；振幅特性在阻帶內是等波紋的，在通帶內是單調下降的切比雪夫II型。(2) 雙線性變換法4工作原理 雙線性變換中數字域頻率和模擬頻率之間的非線性關系限制了它的應用范圍，只有當非線性失真是允許的或能被忽略時，才能采用雙線性變換法，通常低通、高通、帶通和帶阻等濾波器等具有分段恒定的頻率特性，可以采用預畸變的方法來補償頻率畸變，因此可以采用雙線性變換設計方法。(3) 脈沖響應不變法4工作原理 沖激響應不變法遵循的準則是使數字濾波器的單位取樣響應與參照的模擬濾波器的脈沖響應的取樣值完全一樣，即h(n)=ha(nT),其中T為取樣周期。實際是由模擬濾波器轉換成

14、為數字濾波器，就是要建立模擬系統函數Ha（S）與數字系統函數H（z）之間的關系。脈沖響應不變法是從S平面映射到z平面，這種映射不是簡單的代數映射，而是S平面的每一條寬為2/T的橫帶重復地映射到整個z平面。 2.2.2 FIR設計原理由于IIR數字濾波器能夠保留一些模擬濾波器的優良特性，因此應用很廣。但是這些特性是以犧牲線性相位頻率特性為代價的，即用Butterworth、切比雪夫和橢圓法設計的數字濾波器逼近理想的濾波器的幅度頻率特性，得到的濾波器往往是非線性的。在許多電子系統中，對幅度頻率特性和線性相位特性都有較高的要求，所以IIR濾波器在這些系統中往往難以勝任。 設計程序的調試和運行結果3.

15、1聲音信號的提取與加噪試和運行結果3.1.1提取聲音信號x1=wavread(C:UsersacerDesktop鄧紫棋-你把我灌醉.wav,10);x1=x1;x1=x1(1,:);x2=1 2.43 6.17 12.93 22.17 32.25 40.88 45.87 45.87 40.88 32.25 22.17 12.93 6.17 2.43 1;subplot(7,1,1);stem(x1);ylabel(x1(n);title(x1(n);subplot(7,1,2);stem(x2);ylabel(x2(n);title(x2(n);y=conv(x1,x2); subplot(

16、7,1,3);stem(y);ylabel(y);title(x1(n)與x2(n)的卷積);N1=length(x1);N2=length(x2);N=N1+N2-1;X1=fft(x1,N);X2=fft(x2,N);subplot(7,1,4);stem(X1);ylabel(X1);title(x1(n)的N點DFT);subplot(7,1,5);stem(X2);ylabel(X2);title(x2(n)的N點DFT);Y1=X1.*X2;subplot(7,1,6);stem(Y1);ylabel(Y1);title(X1與X2相乘的結果);Y2=ifft(Y1);subplo

17、t(7,1,7);stem(Y2);ylabel(Y2);title(Y1的IDFT結果);運行結果如圖： 圖4 原始信號波形. 對聲音信號加噪MATLAB程序如下：y,fs,bits=wavread(C:UsersacerDesktop鄧紫棋-你把我灌醉.wav);%x：語音數據；fs：采樣頻率；bits：采樣點數 sound(y,fs,bits); %話音回放 n=length (y); %求出語音信號的長度 Y=fft(y,n); %傅里葉變換subplot(2,1,1); plot(y); title(原始信號波形); subplot(2,1,2); plot(abs(Y); titl

18、e(原始信號頻譜); L=length(y);noise=0.5*randn(L,2);y_z=y+noise;sound(y_z,fs,bits);m=length (y_z); %求出被污染語音信號的長度 Y_Z=fft(y_z,n); %傅里葉變換subplot(3,1,1); plot(y_z); title(被污染信號波形); subplot(3,1,2); plot(abs(Y_Z); title(被污染的語音信號頻譜);運行結果如圖： 圖5 被污染的信號3.2 切比雪夫低通濾波器程序的調試和運行結果 如圖所示為原語音信號的時域圖和頻譜圖。在MATLAB中通過wavresd(fil

19、ename)讀取語音信號數據。其程序如下： y,fs,nbits=wavread(C:UsersacerDesktop鄧紫棋-你把我灌醉.wav) 圖6 原始信號 如圖所示為加入噪聲后信號的時域圖和頻譜圖。其主要程序如下：x1=x+zs 圖7 加入噪聲后信號 如圖所示為低通濾波器的幅頻圖。設計指標：wpz=1000Hz; wsz=1200Hz; rp=1; rs=100;其設計程序如下：圖8 切比雪夫低通濾波器如圖所示為濾波后語音信號的時域圖和頻譜圖。其濾波程序如下：yd=filter(az,bz,x1); 圖9 濾波后信號3.3切比雪夫高通濾波器程序的調試和運行結果如圖所示為原語音信號的時域

20、圖和頻譜圖。在MATLAB中通過wavplay()讀取語音信號數據。其程序如下：x,fs,nbits=wavread(C:UsersacerDesktop鄧紫棋-你把我灌醉.wav) 圖10 原語音信號的時域圖和頻譜圖如圖所示為被污染語音信號的時域圖和頻譜圖。其加噪信號程序如下：x1=x+zs1; 圖11 被污染語音信號的時域圖和頻譜圖如圖所示為高通濾波器的幅頻圖。設計指標：wpz=500Hz; wsz=700Hz; rp=1; rs=100;其設計程序如下：wpz=2*700/fs;wsz=2*500/fs;wp=2*fs*tan(wpz*pi/2);ws=2*fs*tan(wsz*pi/2

21、);rp=1;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,az,16384,fs); 圖12 高通濾波器的幅頻圖 如圖所示為濾波后語音信號的時域圖和頻譜圖。其濾波程序如下：yd=filter(az,bz,x1); 圖13 濾波后語音信號的時域圖和頻譜圖 3.4 切比雪夫帶通程序的調試和運行結果如圖所示為原語音信號的時域圖和頻譜圖。在MATLAB中通過wavplay(filename)讀取語音信號數據。其程序如下：x,Fs,bits=wavrea

22、d(C:UsersacerDesktop64ca9e36e190eb780c00b59c4459a6d5.wav); 圖14 原語音信號的時域圖和頻譜圖如圖所示為被污染語音信號的時域圖和頻譜圖。其加噪信號程序如下：x1=x+zs1; 圖15 被污染語音信號的時域圖和頻譜圖 如圖所示為帶通濾波器的幅頻圖。設計指標：fp1=1200 Hz，fp2=3000 Hz，fc1=1000 Hz，fc2=3200 Hz，As=100dB，Ap=1dB;其設計程序如下：wpz=2*1200/fs,2*3000/fs;wsz=2*1000/fs,2*3200/fs;wp=2*fs*tan(wpz*pi/2);w

23、s=2*fs*tan(wsz*pi/2);rp=1;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); 圖16 帶通濾波器的幅頻圖如圖所示為濾波后語音信號的時域圖和頻譜圖。其濾波程序如下：yd=filter(az,bz,x1); 圖17 濾波后語音信號的時域圖和頻譜圖. hamming低通程序的調試和運行結果如圖所示為原語音信號的時域圖和頻譜圖。在MATLAB中通過wavplay(filename)讀取語音信號數據。其程序如下：y,f

24、s,nbits=wavread(C:UsersacerDesktop鄧紫棋-你把我灌醉.wav); 圖18 原語音信號的時域圖和頻譜圖如圖所示為被污染語音信號的時域圖和頻譜圖。其加噪信號程序如下：z=x+y; 圖19 被污染語音信號的時域圖和頻譜圖如圖所示為低通濾波器的幅頻圖和相頻圖。設計指標：wpz=1000Hz; wsz=1200Hz; rs=100;其設計程序如下：wpz=2*1000*pi/fs;wsz=2*1200*pi/fs;rs=100;Bt=wsz-wpz;N0=ceil(6.6*pi/Bt);N=N0+mod(N0+1,2);wc=(wpz+wpz)/2/pi;hn=fir1

25、(N-1,wc,hamming(N);figure(3)freqz(hn,1,32678); 圖20 低通濾波器的幅頻圖和相頻圖如圖所示為濾波后語音信號的時域圖和頻譜圖。其濾波程序如下：hn=filter(hn,1,z); 圖21 濾波后語音信號的時域圖和頻譜圖. hamming高通程序的調試和運行結果如圖所示為原語音信號的時域圖和頻譜圖。在MATLAB中通過wavplay(filename)讀取語音信號數據。其程序如下：y,fs,nbits=wavread(C:UsersacerDesktop鄧紫棋-你把我灌醉.wav); 圖22 原語音信號的時域圖和頻譜圖如圖所示為被污染語音信號的時域圖和

26、頻譜圖。其加噪信號程序如下：z=x+y; 圖23 被污染語音信號的時域圖和頻譜圖如圖所示為高通濾波器的幅頻圖和相頻圖。設計指標：wpzl=2800Hz; wszl=3200Hz; rs=100;其設計程序如下：wpz=2*3000*pi/fs;wsz=2*2800*pi/fs;rs=150;Bt=wpz-wsz;N0=ceil(6.6*pi/Bt);N=N0+mod(N0+1,2);wc=(wpz+wpz)/2/pi;hn=fir1(N-1,wc,high,hamming(N); 圖24 高通濾波器的幅頻圖和相頻圖如圖所示為濾波后語音信號的時域圖和頻譜圖。其濾波程序如下：hn=filter(H

27、n,1,z); 圖25 濾波后語音信號的時域圖和頻譜圖. hamming帶通程序的調試和運行結果如圖所示為原語音信號的時域圖和頻譜圖。在MATLAB中通過wavplay(filename)讀取語音信號數據。其程序如下：y,fs,nbits=wavread(C:UsersacerDesktop鄧紫棋-你把我灌醉.wav); 圖26 原語音信號的時域圖和頻譜圖如圖所示為被污染語音信號的時域圖和頻譜圖。其加噪信號程序如下：t=0:0.000045:0.90855;x1=0.1*sin(2*pi*10*t)+0.2*sin(2*pi*20*t)+0.07*sin(2*pi*6000*t);x=rot9

28、0(x1,2),rot90(x1,2);x=x(1:32768);z=x+y; 圖27 被污染語音信號的時域圖和頻譜圖如圖所示為帶通濾波器的幅頻圖和相頻圖。設計指標：wpz=4500Hz; wsz=4700Hz; wpzl=700Hz; wszl=500Hz; rp=20; rs=50;其設計程序如下：wpzl=1200*pi/fs;wpzu=3000*pi/fs;wszl=1000*pi/fs;wszu=3200*pi/fs;rs=100;Bt=wpzl-wszl;N0=ceil(6.6*pi/Bt);N=N0+mod(N0+1,2);wc=(wpzl+wszl)/2/pi,(wpzu+ws

29、zu)/2/pi;Hn=fir1(N-1,wc,hamming(N); 圖28 帶通濾波器的幅頻圖和相頻圖如圖所示為濾波后語音信號的時域圖和頻譜圖。其濾波程序如下：hn=filter(hn,1,z); 圖29 濾波后語音信號的時域圖和頻譜圖 總結及心得體會近一周的課程設計終于接近尾聲，在這一周里，我感覺到了自己很多的不足，同時也有著諸多的辛酸與收獲。首先，這次的課程設計與以往的課程設計有著很大的不同之處，在時間上相對于以往兩周的課設，這次的時間僅僅為一個星期，或者說4天，這個對我們來說時間上還是比較緊的，需要合理的安排分配時間才能夠及時的完成；在分組上，此次的設計采用一人一組的形式，和以往一組

30、23人的形式相比，這次的實習就缺乏了組員之間的分工與討論，對個人獨立解決問題的能力還是要求相對比較高的。總的來說，這次的實習對我來說是一項挑戰吧，需要我迎難而上認真獨立完成。在課設的過程中讓我感受很深的還是自己能力的不足，讓自己覺得挺失望的。這次的實習主要是以數字信號處理的知識為理論基礎，以MATLAB編程為工具完成對信號的處理。可是自我反省，覺得首先在數字信號的理論基礎這一點上，那些理論知識基本上是學過之后便忘記了，不會付諸于應用，然后再課設的過程中需要自己去花費大量的時間學習，所以就會想要啥自己的理論知識足夠扎實就好了。其次對于MATLAB的應用，覺得以前對于MATLAB的學習還是不夠深入

31、，基本上是依樣畫葫蘆的設計一些較為基礎簡單的程序，不能靈活的運用于自我獨立編程。可是現在的課設和以往的平時的學習就不一樣了，需要自己靈活的運用所學的知識進行編程，不得不說讓我覺得很燒腦，有時候一個程序調試了大半天都沒有得到想要的結果，那種感覺確實很不好。雖說課設燒腦那是肯定的，但是收獲也不少。我覺得課設的目的之一就是以課程設計來實踐完善平時理論知識的不足，同時鍛煉動手操作能力，可以說課程設計和理論學習是一個相輔相成的過程。我也確實對MATLAB這樣一個強大的工具有了更深刻的理解，比如其對聲音信號的提取，加入噪聲的過程以及濾除噪聲的過程。若聯系實際生活我想這與我們生活中的聲音、圖像等方面信息的處

32、理也是有著很大的關聯的，不得不說其實在大學里若認真去學，還是可以學到很多精彩的東西。我很清楚的知道自己是本科畢業以后將要就要面臨找工作的人，通過這樣一次課程設計來補充自己理論知識的不足，提高自己的動手操作能力，為自己以后的就業夯實基礎，我覺得這是一個很有效的途徑。雖然自己現在的能力與以后就業所需要的能力相比還存在著較大的差距，但是我也不可以因此而自暴自棄，不能讓自己的能力不足成為我停滯不前的借口。學習是一個終身的的過程，雖然這次的課設已經結束，我想其在我以后學習生活道路上的影響必定是深遠的。參考文獻1 劉衛國MATLAB程序設計應用M北京：高等教育出版社，20082 趙紅怡數字信號處理及MAT

33、LAB實現M北京：化學工業出版社，20023 高西全，丁玉美數字信號處理（第三版）M陜西：西安電子科技大學出版社，20104 何強，何英MATLAB擴展編程M北京：清華大學出版社，2002 附錄A 切比雪夫濾波器完整程序I切比雪夫低通濾波器x,Fs,bits=wavread(C:UsersacerDesktop鄧紫棋-你把我灌醉.wav);x=x(:,1);figure(1);subplot(2,1,1);plot(x);sound(x,Fs,bits);title(語音信號時域波形圖);y=fft(x,3260);f=(Fs/1630)*1:1630;subplot(2,1,2);plot(

34、f(1:1630),abs(y(1:1630);title(語音信號頻譜圖);t=0:length(x)-1;zs=0.05*cos(2*pi*10000*t/22050);zs0=0.05*cos(2*pi*10000*t/22050000);figure(2);subplot(2,1,1)plot(zs0)title(噪聲信號波形);zs1=fft(zs,1200);sound(zs,Fs,bits); %回放噪音 subplot(2,1,2)plot(f(1:600),abs(zs1(1:600);title(噪聲信號頻譜);%加噪聲之后的語音x1=x+zs; %sound(x1,Fs,

35、bits); 回放加入噪聲后的語音y1=fft(x1,1200);figure(3);subplot(3,1,1);plot(x1);title(加入噪聲后的信號波形);subplot(3,1,2);plot(f(1:600),abs(y1(1:600);title(加入噪聲后的信號頻譜)%低通濾波器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=35;N,wc=cheb2ord(wp,ws,rp,rs,s);b,a=cheby2(N,rs,wc,s);bz,az=bilinear(

36、b,a,fs);H,W=freqz(bz,az);figure(4); plot(w,abs(h);title(切比雪夫低通濾波器);xlabel(頻率（HZ）);ylabel(耗損（dB）);grid on;yd=filter(az,bz,x1);figure(5);subplot(2,1,1);plot(yd);title(濾波后的時域波形圖);ydd=fft(yd,800);subplot(2,1,2);plot(f(1:600),abs(ydd(1:600);title(濾波后的頻域波形圖);sound(yd,Fs,bits);II切比雪夫高通濾波器x,Fs,bits=wavread(

37、C:UsersacerDesktop鄧紫棋-你把我灌醉.wav);x=x(:,1);figure(1);subplot(2,1,1);plot(x);sound(x,Fs,bits);title(語音信號時域波形圖);y=fft(x,3260);f=(Fs/1630)*1:1630;subplot(2,1,2);plot(f(1:1630),abs(y(1:1630);title(語音信號頻譜圖);%被污染的語音L=length(x);zs0=0.5*randn(L,2);figure(2);subplot(2,1,1)plot(zs0)title(噪聲信號波形);zs1=fft(zs,120

38、0);sound(zs,Fs,bits); %回放噪音 subplot(2,1,2)plot(f(1:600),abs(zs1(1:600);title(噪聲信號頻譜);%加噪聲之后的語音x1=x+zs;%sound(x1,Fs,bits); %回放加入噪聲后的語音y1=fft(x1,1200);figure(3);subplot(3,1,1);plot(x1);title(加入噪聲后的信號波形);subplot(2,1,2);plot(f(1:600),abs(y1(1:600);title(加入噪聲后的信號頻譜);%高通fp=3000;fs=2800;Fs=22050;rp=1;rs=10

39、0;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs=Fs/Fs;wap=2*tan(wp/2);was=2*tan(ws/2);N,wc=cheb1ord(wp,ws,rp,rs,s);B,A=cheby1(N,rp,wc,high,s);Bz,Az=bilinear(B,A,Fs);figure(4);h,w=freqz(Bz,Az,512,Fs*22050);plot(w,abs(h);title(切比雪夫高通濾波器);xlabel(頻率（HZ）);ylabel(耗損（dB）);grid on;yd=filter(Bz,Az,x1);figure(5);subplot(3,1

40、,2);plot(yd);title(濾波后的時域波形圖);ydd=fft(yd,800);subplot(3,1,3);plot(f(1:600),abs(ydd(1:600);title(濾波后的頻域波形圖);切比雪夫帶通濾波器x,Fs,bits=wavread(C:UsersacerDesktop64ca9e36e190eb780c00b59c4459a6d5.wav);x=x(:,1);figure(1);subplot(2,1,1);plot(x);sound(x,Fs,bits);title(語音信號時域波形圖);y=fft(x,3260);f=(Fs/1630)*1:1630;s

41、ubplot(2,1,2);plot(f(1:1630),abs(y(1:1630);title(語音信號頻譜圖);%被污染的語音t=0:length(x)-1;zs=0.05*cos(2*pi*10000*t/22050);zs0=0.05*cos(2*pi*10000*t/22050000);figure(2);subplot(2,1,1)plot(zs0)title(噪聲信號波形);zs1=fft(zs,1200);sound(zs,Fs,bits); %回放噪音 subplot(2,1,2)plot(f(1:600),abs(zs1(1:600);title(噪聲信號頻譜);%加噪聲之

42、后的語音x1=x+zs;%sound(x1,Fs,bits); %回放加入噪聲后的語音y1=fft(x1,1200);figure(3);subplot(2,1,1);plot(x1);title(加入噪聲后的信號波形);subplot(2,1,2);plot(f(1:600),abs(y1(1:600);title(加入噪聲后的信號頻譜);fs1=500;fp1=800;fpu=1500;fsu=2000;Fs=22050rp=1;rs=100;wp=2*pi*fp1,fpu/Fs;ws=2*pi*fs1,fsu/Fs;Fs=Fs/Fs;wap=2*tan(wp/2);was=2*tan(w

43、s/2);N,wc=cheb1ord(wp,ws,rp,rs,s);B,A=cheby1(N,rp,wc,s);Bz,Az=bilinear(B,A,Fs);h,w=freqz(Bz,Az,512,Fs*22050);figure(4);plot(w,abs(h);title(切比雪夫帶通濾波器);xlabel(頻率（HZ）);ylabel(耗損（dB）);grid on;yd=filter(Bz,Az,x1);figure(5);subplot(2,1,1);plot(yd);title(濾波后的時域波形圖);ydd=fft(yd,800);subplot(2,1,2);plot(f(1:6

44、00),abs(ydd(1:600);title(濾波后的頻域波形圖);sound(yd,Fs,bits);附錄加窗濾波器完整程序I hamming窗低通濾波器完整程序fs=22050;y,fs,nbits=wavread(C:UsersacerDesktop鄧紫棋-你把我灌醉.wav);y=y(1:32768)+0.5;figure(1)subplot(2,1,1);plot(abs(y);xlabel(n);ylabel(y);title(原語音信號時域波形)y=y-0.5;Y=fft(y,32768);f=22050*(0:16384)/32768;subplot(2,1,2);plot

45、(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,2),rot90(x1,2);x=x(1:32768);z=x+y;figure(2)subplot(2,1,1);plot(z);xlabel(n);ylabel(z);title(被污染語音信號時域波形)Z=fft(z,32768);f=22050*(0:16384)/32768;subplot(2,1,2);plot(f,abs

46、(Z(1:16385);xlabel(f);ylabel(Z);title(被污染語音信號頻譜圖)wpz=2*1000*pi/fs;wsz=2*1200*pi/fs;rs=100;Bt=wsz-wpz;N0=ceil(6.6*pi/Bt);N=N0+mod(N0+1,2);wc=(wpz+wpz)/2/pi;hn=fir1(N-1,wc,hamming(N);figure(3)freqz(hn,1,32678);figure(4)hn=filter(hn,1,z);H=fft(hn,32768);f=22050*(0:16382)/32768;sound(H,fs,nbits);subplot

47、(2,1,2);plot(f,abs(H(1:16383);xlabel(f);ylabel(H);title(濾波后語音信號頻譜圖)subplot(2,1,1);hn=hn+0.5;plot(abs(hn);xlabel(n);ylabel(hn);title(濾波后語音信號時域波形)hamming高通濾波器完整程序fs=22050;y,fs,nbits=wavread(C:UsersacerDesktop鄧紫棋-你把我灌醉.wav);y=y(1:32768);figure(1)subplot(2,1,1);plot(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.1*sin(2*pi*100*t)+0.2*sin(2*pi*220*t);x=rot90(x1,2),rot90(x1,