語音信號的數(shù)字濾波處理(五)

語音信號的數(shù)字濾波處理圖1所示，程序如下：n=0:50;A=10;a=10*sqrt(2)*pi;w0=10*sqrt(2)*pi;t0=0.001;t1=0.020;t2=0.005;x=A*exp(-a*n*t0).*sin(w0*n*t0);y1=A*exp(-a*n*t1).*sin(w0*n*t1);y2=A*exp(-a*n*t2).*sin(w0*n*t2);subplot(321);stem(n,x,'r'),gridontitle('離散時間信號')xlabel('n');ylabel('幅度')subplot(322);plot(n,x,'g'),gridontitle('連續(xù)時間信號')xlabel('t');ylabel('幅度')subplot(323);stem(n,y1),gridontitle('50Hz連續(xù)時間信號')xlabel('n');ylabel('幅度')subplot(324);stem(n,y2,'m'),gridontitle('200Hz理想采樣信號序列');xlabel('n');ylabel('幅度')subplot(325);stem(n,x,'k'),gridontitle('1000Hz連續(xù)時間信號')xlabel('n');ylabel('幅度')t=0:0.001:1;xx=A*exp(-a*n*t0).*sin(w0*n*t0);X=fft(xx);subplot(326)plot(abs(X));xlabel('f/Hz');ylabel('|X(ejw)|')gridon；圖2.4采樣定理圖示3課程設(shè)計的基本原理3.1IIR濾波器的設(shè)計原理IIR數(shù)字濾波器的設(shè)計一般是利用目前已經(jīng)很成熟的模擬濾波器的設(shè)計方法來進行設(shè)計，通常采用模擬濾波器原型有巴特沃斯濾波器、切比雪夫濾波器、橢圓濾波器等。IIR數(shù)字濾波器單位取樣響應(yīng)是時寬無限的h(n)，即n0＜n＜∞。數(shù)字濾波器按照實現(xiàn)的方法和結(jié)構(gòu)形式分為遞歸型或非遞歸型兩類。遞歸型數(shù)字濾波器的當(dāng)前輸出y(n)是輸入x(n)的當(dāng)前值和以前各輸入值x(n)，x(n–1)，….，及以前各輸出值y(n)，y(n–1)，….的函數(shù)。一個N階遞歸型數(shù)字濾波器(IIR濾波器)的差分方程如公式2-1所示：2-1其中，式中的系數(shù)至少有一項不為零。說明必須將延時的輸出序列進行反饋。遞歸系統(tǒng)的傳統(tǒng)函數(shù)定義如公式2-2所示：2-2遞歸系統(tǒng)的傳遞函數(shù)H(z)在Z平面上不僅有零點，而且有極點。非遞歸型數(shù)字濾波器當(dāng)前的輸出值y(n)僅為當(dāng)前及以前的輸入序列的函數(shù)，而與以前的各個輸出值無關(guān)，因此從結(jié)構(gòu)上看非遞歸系統(tǒng)沒有反饋環(huán)路。3.2FIR濾波器的設(shè)計原理FIR數(shù)字濾波器的單位取樣響應(yīng)是時寬有限的h(n)，即n1≤n＜n2。一個N階的非遞歸型數(shù)字濾波器(FIR濾波器)的差分方程如公式2-3所示：2-3差分方程式中的系數(shù)等于單位取樣響應(yīng)的序列值h(n)，其系統(tǒng)函數(shù)H(z)可以表示為公式2-4所示的形式：2-4H(z)是的多項式，因此它的極點只能在Z平面的原點上。FIR濾波器通常采用窗函數(shù)方法來設(shè)計。窗設(shè)計的基本思想是，首先選擇一個適當(dāng)?shù)睦硐脒x頻濾波器（它總是具有一個非因果，無限持續(xù)時間脈沖響應(yīng)），然后街區(qū)（加窗）它的脈沖響應(yīng)得到線性相位和因果FIR濾波器。我們用表示理想的選頻濾波器，它在通帶上具有單位增益和線性相位，在阻帶上具有零響應(yīng)。一個帶寬wc<pi的低通濾波器由下式給定：為了從hd(n)得到一個FIR濾波器，必須同時在兩邊截取hd（n）。而要得到一個因果的線性相位濾波器，它的h(n)長度為N，必須滿足公式2-5：2-5h(n)可以看做是hd(n)與窗函數(shù)w（n）的乘積，如公式2-6所示：2-6這種操作叫做加窗，其中根據(jù)w（n）的不同定義，可以得到不同的窗結(jié)構(gòu)。在頻域中，因果FIR濾波器響應(yīng)由和窗響應(yīng)的周期卷積得到，如公式2-7所示:2-7常用的窗函數(shù)有矩形窗、漢寧（HANNING）窗、海明（矩形）窗、布萊克曼（BLACKMAN）窗、凱澤（KAISER）窗等。4設(shè)計程序的調(diào)試及運行結(jié)果4.1IIR濾波器的設(shè)計(1)選擇來自于D盤的“ding.wav”聲音作為語音信號(2)根據(jù)不同的濾波器選擇不同的噪聲信號加到語音信號中，得到被污染的語音信號。(3)分別設(shè)計巴特沃思濾波器（低通、帶通、高通）以及Kaiser窗濾波器（低通、帶通、高通）對被污染的語音信號濾波，濾掉相應(yīng)的噪音信號，得到符合要求的語音信號。方案設(shè)計流程圖如圖4.1所示。（圖為巴特沃思濾波器，Kaiser窗濾波器與巴特沃思類似）4.1總體方案設(shè)計流圖4.11巴特沃斯低通濾波器的設(shè)計(1)添加高頻的噪聲信號zs=0.05*cos(2*pi*10000*t/22050)，MATLAB程序如下。噪聲信號的時域及幅頻圖如圖3.2所示。%噪聲信號zt=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('噪聲信號頻譜');圖4.2噪聲信號時域和頻譜圖(2)設(shè)計巴特沃思低通濾波器，MATLAB程序如下。濾波器圖如圖3.3所示。%低通濾波fp=3000;fs=3500;Fs=22050;rp=1;rs=10;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs1);figure(4);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯低通濾波器');xlabel('頻率（HZ）');ylabel('耗損（dB）');gridon;圖4.3巴特沃思低通濾波器圖4.12巴特沃斯高通濾波器的設(shè)計(1)添加低頻噪聲信號zs0=0.05*cos(2*pi*100*t/22050)，MATLAB程序如下。噪聲信號的時域及幅頻圖如圖3.4所示。t=0:length(x)-1;zs0=0.05*cos(2*pi*100*t/22050);figure(2);subplot(2,1,1)plot(zs0)title('噪聲信號波形');zs1=fft(zs0,1200);%sound(zs,FS,bits);%回放噪音subplot(2,1,2)plot(f(1:600),abs(zs1(1:600)));title('噪聲信號頻譜');圖4.4噪聲信號時域和頻譜圖(2)設(shè)計巴特沃思高通濾波器，MATLAB程序如下。濾波器圖如圖3.5所示。%高通濾波fp=600;fs=400;Fs=22050;rp=1;rs=10;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'high','s');[Bz,Az]=bilinear(B,A,Fs1);figure(4);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯高通濾波器');xlabel('頻率（HZ）');ylabel('耗損（dB）');gridon;圖4.5巴特沃思高通濾波器圖4.13巴特沃斯高通濾波器的設(shè)計(1)添加低頻噪聲信號zs1=0.05*cos(2*pi*100*t/22050)，MATLAB程序如下。噪聲信號的時域及幅頻圖如圖3.6所示。t=0:length(x)-1;zs0=0.05*cos(2*pi*100*t/22050);figure(2);subplot(2,1,1)plot(zs0)title('噪聲信號波形');zs1=fft(zs0,1200);%sound(zs,FS,bits);%回放噪音subplot(2,1,2)plot(f(1:600),abs(zs1(1:600)));title('噪聲信號頻譜');圖4.6噪聲信號時域和頻譜圖(2)設(shè)計巴特沃思帶通濾波器，MATLAB程序如下。濾波器圖如圖3.7所示。%帶通濾波fp=[600,6000];fs=[400,7000];Fs=22050;rp=1;rs=10;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs1);figure(4);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯帶通濾波器');xlabel('頻率（HZ）');ylabel('耗損（dB）');gridon;圖4.7巴特沃思帶通濾波器圖4.2FIR數(shù)字濾波器的設(shè)計4.21Kaiser窗低通程序的設(shè)計技術(shù)指標：fr=1000。源程序如下。仿真波形見圖4.21。[y,fs,bit]=wavread('1');fr=1000;wr=2*pi*fr;Window=kaiser(8,beta);b=fir1(7,fr/(fs/2),Window);[h,w]=freqz(b,1);subplot(221);plot(w*fs/(2*pi),abs(h));title('低通濾波');y1=filter(b,1,y);y2=fftfilt(b,y);subplot(222);plot(y);title('原始語音信號');subplot(223);plot(y1);title('IIR濾波后語音信號');subplot(224);plot(y2);title('FIR濾波后語音信號');sound(y1,fs,bit);圖4.21凱塞窗低通濾波器仿真圖4.22Kaiser窗帶通濾波器技術(shù)指標：fr=1000,fp=3000。源程序如下。仿真波形見圖4.22。[y,fs,bit]=wavread('1');fr=1000;fp=3000wr=2*pi*fr;Window=kaiser(16,beta);b=fir1(15,[fr/(fs/2),fp/(fs/2)],Window);[h,w]=freqz(b,1);subplot(221);plot(w*fs/(2*pi),abs(h));title('帶通濾波');y1=filter(b,1,y);y2=fftfilt(b,y);subplot(222);plot(y);title('原始語音信號');subplot(223);plot(y1);title('IIR濾波后語音信號');subplot(224);plot(y2);title('FIR濾波后語音信號');sound(y1,fs,bit);圖4.22凱塞窗帶通濾波器仿真波形4.23Kaiser窗高通濾波器指標：fr=3000。源程序如下。仿真波形見圖4.23。[y,fs,bit]=wavread('1');fr=3000;wr=2*pi*fr;Window=kaiser(8,beta);b=fir1(7,fr/(fs/2),'high');[h,w]=freqz(b,1);subplot(221);plot(w*fs/(2*pi),abs(h));title('低通濾波');y1=filter(b,1,y);y2=fftfilt(b,y);subplot(222);plot(y);title('原始語音信號');subplot(223);plot(y1);title('IIR濾波后語音信號');subplot(224);plot(y2);title('FIR濾波后語音信號');sound(y1,fs,bit);圖4.23凱塞窗高通濾波器仿真波形4.3回放語音信號經(jīng)過以上的加噪聲處理后，可在Matlab中用函數(shù)sound對聲音進行回放。其調(diào)用格式：sound(y,Fs)，sound(y)和sound(y,Fs,bits)。可以察覺濾波前后的聲音有明顯的變化。5心得與體會在本次課程設(shè)計的過程中，主要遇到的問題是原始音樂信號的選取。有些音樂信號在MATLAB中運行時出錯，開始以為自己的操作問題，在跟同學(xué)商討多次更換音樂信號后得到了解決。其次就是在濾波器參數(shù)的選取上，開始時不清楚怎么選取，又在課本中找了好久才找到。學(xué)習(xí)的過程是相互討論共同進步的，多多討論課題中遇到的問題，可以鞏固我們的知識掌握能力，增加熟練運用度。從本次課程設(shè)計的中心來看，課題是希望將數(shù)字信號處理技術(shù)應(yīng)用于某一實際