《數(shù)字信號(hào)處理》課程研究性學(xué)習(xí)報(bào)告_DSP3

1、數(shù)字信號(hào)處理課程研究性學(xué)習(xí)報(bào)告姓名 張卓 學(xué)號(hào) 09271119 同組成員 王欣然 王一明 馬高飛 郭鵬飛 指導(dǎo)教師 胡健 時(shí)間 2011年5月29日 數(shù)字信號(hào)處理綜合應(yīng)用專題研討【目的】(1) 能夠靈活應(yīng)用DFT（FFT）分析實(shí)際信號(hào)的頻譜。(2) 熟悉通過IIR和FIR數(shù)字濾波器進(jìn)行實(shí)際系統(tǒng)設(shè)計(jì)的方法。(3) 培養(yǎng)學(xué)生自主學(xué)習(xí)能力，以及發(fā)現(xiàn)問題、分析問題和解決問題的能力。【研討題目】 1附件給出了一段含有噪聲的音頻信號(hào)。（1）分析該信號(hào)的頻譜特點(diǎn)。（2）通過脈沖響應(yīng)不變法設(shè)計(jì)一個(gè)濾波器對(duì)其進(jìn)行處理，得到有用信息，自主確定各項(xiàng)指標(biāo)。（3）使用雙線性變換法重做（2）。（4）試定量比較上述兩種濾

2、波器的各項(xiàng)性能，畫出能說明性能差異的相關(guān)圖形，對(duì)比并解釋。（5）通過不同窗函數(shù)法設(shè)計(jì)FIR濾波器對(duì)其進(jìn)行處理，比較不同窗函數(shù)的濾波效果，給出分析和結(jié)論。 (6) 請(qǐng)嘗試采用其它的音頻信號(hào)，混入不同的噪聲，利用所學(xué)的濾波方法進(jìn)行分析，會(huì)得到什么樣的效果？【題目分析】本題討論用IIR和FIR數(shù)字濾波器進(jìn)行實(shí)際系統(tǒng)設(shè)計(jì)的方法?！痉抡娼Y(jié)果】(1)(2)脈沖響應(yīng)不變法雙線性變換法(3)Hamming窗Blackman窗Kaiser窗(4)Butterworth(5)【結(jié)果分析】對(duì)幾種濾波器得到的音頻信號(hào)進(jìn)行對(duì)比。原始音頻信號(hào)頻譜在頻率為104Hz以及4*104Hz處分別有兩處毛刺，關(guān)于2.6*104Hz

3、大致對(duì)稱。對(duì)于IIR濾波器，用Butterworth濾波器，指標(biāo)：fp=6000Hz，fs=8000Hz，Ap=0.5dB，As=50dB。對(duì)于FIR濾波器，用Hamming窗，指標(biāo)：fp=6000;fs=8000。用種濾波器濾波后均去處了毛刺，都達(dá)到了濾波的效果。畫出兩種濾波器的幅度相應(yīng)圖形，可以看出，Butterworth濾波器的幅度相應(yīng)較為平滑，沒有等波紋振蕩，通帶比FIR濾波器的小，過渡帶也比FIR濾波器要小。而FIR濾波器的幅度響應(yīng)在阻帶有所波動(dòng)。加Hamming窗的濾波器N=79大于Butterworth濾波器的N,IIR數(shù)字濾波器可以用較少的階數(shù)獲得較好的幅度相應(yīng)，但由于其結(jié)構(gòu)存

4、在反饋，在實(shí)現(xiàn)時(shí)可能會(huì)造成系統(tǒng)不穩(wěn)定，并且是非線性相位的。由于FIR數(shù)字濾波器的單位脈沖相應(yīng)是有限長(zhǎng)，因而系統(tǒng)是穩(wěn)定的，并且很容易設(shè)計(jì)成線性相位?！咀灾鲗W(xué)習(xí)內(nèi)容】一、IIR1. 模擬模擬低通濾波器（1）Butterworth型低通濾波器N,wc=buttord(wp,ws,Ap,As,s)num,den=butterN,wc,s（2）Chebyshev I型低通濾波器N,wc=cheb1ord(wp,ws,Ap,As,s)num,den=cheby1N,Ap,wc,s（3）Chebyshev II型低通濾波器N,wc=cheb2ord(wp,ws,Ap,As,s)num,den=cheby2N

5、,As,wc,s（4）橢圓型低通濾波器N,wc=ellipord(wp,ws,Ap,As,s)num,den=ellipN,Ap,As,wc,s2. 模擬頻率變換（1）低通-高通 numt,dent=lp2hp(num,den,w0)（2）低通-帶通 numt,dent=lp2bp(num,den,w0,B)（3）低通-帶阻 numt,dent=lp2bs(num,den,w0,B)3. 脈沖響應(yīng)不變法 numd,dend=impinvar(num,den,Fs)雙線性變換法numd,dend=bilinear(num,den,Fs)二、FIR 窗函數(shù)法w=hamming(N)w=blackm

6、an(N) w=kaiser(N,beta) M,Wc,beta,ftype=kaiserord(f,a,dev) 直接實(shí)現(xiàn) h=firl(M,Wc,ftype,window) h=firl(M,f,a,npt,window)三、給音頻加噪noise=0.04*rand(N,2);【閱讀文獻(xiàn)】1 陳后金，薛健，胡健. 數(shù)字信號(hào)處理M. 第二版. 北京：高等教育出版社. 2008. 146-152，187-193.【發(fā)現(xiàn)問題】 (專題研討或相關(guān)知識(shí)點(diǎn)學(xué)習(xí)中發(fā)現(xiàn)的問題)：我們本來(lái)打算用matlab來(lái)觀察wc的變化是如何對(duì)濾波器的性能造成影響的，在觀察程序的時(shí)候，我們?cè)噲D改變wc的值，但是num,d

7、en=butter(N,wc,'s');numd,dend=bilinear(num,den,f);是一個(gè)固定結(jié)構(gòu)這里返回的參數(shù)N是BW的階數(shù)，wc是3dB截頻，而且在書上我們查閱到wc是由阻帶確定的，所以在現(xiàn)在所查閱資料的情況下我們無(wú)法對(duì)wc進(jìn)行改變。但是我們對(duì)通帶截頻和阻帶截頻進(jìn)行了一定的研究，就是改變兩個(gè)值來(lái)觀察其對(duì)濾波性能的影響。 原始信號(hào) fp=2000,fs=8000 fp=1000,fs=8000 fp=1000,fs=9000從圖中我們發(fā)現(xiàn)，這幾個(gè)濾波器的濾波效果近乎完全相同的，所以裕量只是讓濾波器的性能個(gè)更加完善，但是并不從本質(zhì)上改變?yōu)V波性能！在做FIR濾波器

8、濾波時(shí)，我們經(jīng)過多次嘗試發(fā)現(xiàn)無(wú)論用何種指標(biāo)，濾波后都不能播出音頻，這個(gè)問題一直不能得到解決?！締栴}探究】在IIR濾波器設(shè)計(jì)過程中，由于利用指標(biāo)數(shù)據(jù)的不同，造成裕量的出現(xiàn)。討論利用不同指標(biāo)出現(xiàn)裕量對(duì)濾波器性能的影響，以及如何有效地利用它？在FIR濾波器設(shè)計(jì)過程中，由于所選用窗函數(shù)的不同，導(dǎo)致對(duì)信號(hào)濾波的效果不同。本題的語(yǔ)音信號(hào)加入的是雙頻噪聲，討論如果對(duì)信號(hào)加入其它形式的噪聲，采用何種濾波器更合適。在設(shè)計(jì)過程中，我們選用的是BW濾波器。在BW設(shè)計(jì)過程中我們需要的一個(gè)最重要的參數(shù)是wc，但是我們做計(jì)算的wc是在一個(gè)范圍內(nèi)的，這個(gè)范圍與通帶截頻，阻帶截頻，通帶衰減與阻帶衰減有直接的關(guān)系。wc如果更加

9、偏向于wp，則阻帶就會(huì)裕量較大，但是如果wc更加偏向于ws，則其通帶衰減較大。但是滿足此范圍的wc的都能實(shí)現(xiàn)該濾波器的性能。本題應(yīng)用了三種方法設(shè)計(jì)FIR數(shù)字濾波器，通過三種不同的函數(shù)的求解過程對(duì)初始的噪音信號(hào)進(jìn)行濾波。但由于基本原理都是應(yīng)用FIR加窗的方法進(jìn)行設(shè)計(jì)，我們先后分析了hamming窗、blackman窗、Kaiser窗對(duì)于此例題中的噪音信號(hào)的影響，發(fā)現(xiàn)在參數(shù)Wp、Ws設(shè)定相同的情況下，濾除噪音信號(hào)的效果及產(chǎn)生的頻譜幾乎相同，差別不大?！痉抡娉绦颉浚?）y,f=wavread('nihao'); wavplay(y,f); n = length (y) ; m=fft

10、(y,n);subplot(2,1,1);plot(y);title(' Spectrum of the original signal in time domain ');subplot(2,1,2);plot(abs(m);xlabel('Hz');title(' Spectrum of the original signal in frequency domain ');（2）%默認(rèn)抽樣頻率f=44100 由圖可設(shè)fp=6000Hz fs=8000HzAp=0.5As=50【脈沖響應(yīng)不變法設(shè)計(jì)濾波器】y,f=wavread('nih

11、ao'); wavplay(y,f); n=length(y) ; m=fft(y,n);fp=6000;fs=8000;Ap=0.5;As=50;wp=2*pi*fp;ws=2*pi*fs;omegap=wp/f;omegas=ws/f;N,wc=buttord(wp,ws,Ap,As,'s'); num,den=butter(N,wc,'s');numd,dend=impinvar(num,den,f);z1=filter(numd,dend,y);wavplay(z1,f);m1=fft(z1); subplot(2,1,1);plot(z1);t

12、itle('Spectrum of the original signal in time domain after be filtered');subplot(2,1,2);plot(abs(m1);title('Spectrum of the original signal in frequency domain after be filtered ');【雙線性變換法設(shè)計(jì)的濾波器】y,f=wavread('nihao'); wavplay(y,f); n=length(y) ; m=fft(y,n);fp=6000;fs=8000;Ap=0

13、.5;As=50;wp=2*pi*fp;ws=2*pi*fs;omegap=wp/f;omegas=ws/f;N,wc=buttord(wp,ws,Ap,As,'s'); num,den=butter(N,wc,'s');numd,dend=bilinear(num,den,f);z1=filter(numd,dend,y);wavplay(z1,f);m1=fft(z1); subplot(2,1,1);plot(z1);title('Spectrum of the original signal in time domain after be fil

14、tered');subplot(2,1,2);plot(abs(m1);title('Spectrum of the original signal in frequency domain after be filtered ');（3）【Hamming窗】y,f=wavread('nihao.wav'); wavplay(y,f); n=length(y) ; Y=fft(y,n);fp=6000;fs=8000;Ap=0.5;As=50;wp=2*pi*fp/f;ws=2*pi*fs/f; N=ceil(7*pi/(ws-wp); %求濾波器長(zhǎng)度NN

15、=mod(N+1,2)+N; %為實(shí)現(xiàn)FIR類型I偶對(duì)稱濾波器，應(yīng)確保N為奇數(shù)M=N-1;fprintf('N=%.0fn',N);w=hamming(N)' %使用哈明窗，并將列向量變?yōu)樾邢蛄縒c=(wp+ws)/2; %截止頻率取通阻帶頻率的平均值k=0:M;hd=(Wc/pi)*sinc(Wc*(k-0.5*M)/pi); %建立理想低通濾波器h=hd.*w; %求FIR系統(tǒng)函數(shù)系數(shù)z1=filter(hd,1,y);wavplay(z1,f);magdb=20*log10(abs(z1);subplot(2,1,1);plot(z1);title('Sp

16、ectrum of the signal with Hamming in time domain ');z2=fft(z1,n);subplot(2,1,2);plot(abs(z2);title('Spectrum of the signal with Hamming in freq domain ');【Blackman窗】y,f=wavread('nihao.wav'); wavplay(y,f); n=length(y) ; Y=fft(y,n);fp=6000;fs=8000;Ap=0.5;As=50;wp=2*pi*fp/f;ws=2*pi*

17、fs/f; N=ceil(7*pi/(ws-wp);N=mod(N+1,2)+N;M=N-1;fprintf('N=%.0fn',N);w=blackman(N)'Wc=(wp+ws)/2;k=0:M;hd=(Wc/pi)*sinc(Wc*(k-0.5*M)/pi);h=hd.*w;z1=filter(hd,1,y);wavplay(z1,f);magdb=20*log10(abs(z1);subplot(2,1,1);plot(z1);title('Spectrum of the signal with Blackman in time domain '

18、;);z2=fft(z1,n);subplot(2,1,2);plot(abs(z2);title('Spectrum of the signal with Blackman in freq domain ');ylabel('dB');【Kaiser窗】clc;clear all;y,f=wavread('nihao.wav'); wavplay(y,f); n=length(y) ; Y=fft(y,n);fp=6000*2/f;fs=8000*2/f;Ap=0.5;As=50;Rp=1-10.(-0.05*Ap);Rs=10.(-0.05*

19、As);f1=fp,fs;a=1,0;dev=Rp,Rs;M,Wc,beta,ftype = kaiserord(f1,a,dev);M=mod(M,2)+M;h = fir1(M,Wc,ftype,kaiser(M+1,beta);z1=filter(h,1,y);wavplay(z1,f); magdb=20*log10(abs(z1);subplot(2,1,1);plot(z1);title('Spectrum of the signal with Kaiser in time domain ');z2=fft(z1,n);subplot(2,1,2);plot(abs

20、(z2);title('Spectrum of the signal with Kaiser in freq domain ');（4）%Butterworthy,f=wavread('nihao'); fp=6000;fs=8000;Ap=0.5;As=50;wp=fp*2*pi;ws=fs*2*pi;omegap=wp/f;omegas=ws/f;N,wc=buttord(wp,ws,Ap,As,'s');num,den=butter(N,wc,'s');numd,dend=bilinear(num,den,f);H,w=fr

21、eqz(numd,dend);dbH=20*log(abs(H)/max(abs(H); plot(w/pi,abs(H);%Hamming窗y,f=wavread('nihao'); fp=6000;fs=8000;Ap=0.5;As=50;wp=2*pi*fp/f;ws=2*pi*fs/f; N=ceil(7*pi/(ws-wp);%求濾波器長(zhǎng)度NN=mod(N+1,2)+N;%為實(shí)現(xiàn)FIR類型I偶對(duì)稱濾波器，應(yīng)確保N為奇數(shù)M=N-1;fprintf('N=%.0fn',N);w=hamming(N)'%使用哈明窗，并將列向量變?yōu)樾邢蛄縒c=(wp+

22、ws)/2;%截止頻率取通阻帶頻率的平均值k=0:M;hd=(Wc/pi)*sinc(Wc*(k-0.5*M)/pi);%建立理想低通濾波器h=hd.*w;%求FIR系統(tǒng)函數(shù)系數(shù)omega=linspace(0,pi,512);mag=freqz(h,1,omega);magdb=20*log10(abs(mag);plot(omega/pi,magdb);ylabel('dB');title('hamming');（5）z,f=wavread('paopao');wavplay(z,f);N=length(z);noise=0.04*rand(

23、N,2);X=fft(noise);subplot(2,1,1);plot(noise);title(' Spectrum of the noise signal in time domain ');subplot(2,1,2);plot(abs(X);title(' Spectrum of the noise signal in frequency domain ');axis(0,1200,0,5);figure;s=z+noise;X1=fft(s);wavplay(s,f);subplot(2,1,1);plot(s);title(' Spect

24、rum of the original signal in time domain ');subplot(2,1,2);plot(abs(X1);title(' Spectrum of the original signal in frequency domain ');figure;fp=150;fs=250;Ap=3;As=20;Fs=800;T=1/Fs;wp=fp/Fs*2*pi;ws=fs/Fs*2*pi;Omgp=(2/T)*tan(wp/2);Omgs=(2/T)*tan(ws/2);n,Omgc=buttord(Omgp,Omgs,Ap,As,'s');ba,aa=butter(n,Omgc,'s');bd,ad=bilinear(ba,aa,Fs);z1=filter(bd,ad,s);wavplay(z1,f);m=fft(z1);subplot(2,1,1);plot(z);title('Spectrum of the original s