【文章內(nèi)容簡(jiǎn)介】 )))*225。 figure。 imshow(A)。 在 MATLAB 中執(zhí)行了 FFT 后，使用了 fftshift 函數(shù)調(diào)整，以使頻譜圖像的原點(diǎn)從起始點(diǎn) (0， 0)，移到圖像的中心點(diǎn) (M/2,N/2)，對(duì)應(yīng)的圖 的傅立葉頻譜能量圖。從圖中可以看出：豎直方向理想單一頻率的正弦分布紋理的頻譜能量 集中在水平方向的三個(gè)點(diǎn)。 圖 紋理頻譜圖 (3)離散信號(hào) /序列 離散時(shí)間周期信號(hào)能夠用具有諧波關(guān)系的復(fù)指數(shù)序列的線性組合來(lái)表示，稱(chēng)為離散傅里葉級(jí)數(shù)。將這一概念推廣應(yīng)用到離散時(shí)間非周期信號(hào)，認(rèn)為離散時(shí)間非周期信號(hào)也能夠用具有諧波關(guān)系的復(fù)指數(shù)序列的線性組合來(lái)表示。 當(dāng)離散時(shí)間周期信號(hào)的周期 N 趨于無(wú)窮大時(shí)，則離散時(shí)間周期信號(hào)就轉(zhuǎn)化為離散時(shí)間非周期信號(hào)，其離散頻譜就轉(zhuǎn)化為連續(xù)頻譜，稱(chēng)為離散時(shí)間傅里葉變換(Discrete Time Fourier Transform， DTFT)。 13 以 ( ) 1 2 sin ( 2 * 1 0 ) 3 sin ( 2 * 1 5 )x t t t??? ? ?，采樣頻率 100Hz，采樣 100 個(gè)點(diǎn)，形成離散信號(hào)，利用 MATLAB 進(jìn)行頻譜分析。 程序如下： N=100。Fs=100。t=(0:(N1))/Fs。 xn=1+2*sin(2*pi*10*t)+3*sin(2*pi*15*t)。 XK=abs(fft(xn,N))。f=(0:N/2)*Fs/N。 XK(1)=XK(1)/N。XK(2:(N/2+1))=XK(2:(N/2+1))*2/N。 stem(f,XK(1:(N/2+1)))。axis([1 N/2 0 5])。 grid on。xlabel(39。f(Hz)39。)。ylabel(39。|X( f )|39。)。 程序運(yùn)行結(jié)果如圖 34所示。 圖 N=100 離散信號(hào)頻譜 從圖 可以看出， DFT 法分析的結(jié)果和實(shí)際信號(hào)吻合得很好，說(shuō)明該方法確實(shí)有效。 但是，需要注意的是，離散傅里葉變換在頻域是離散的，即限制在基頻整數(shù)倍上，只能得到信號(hào) K 次諧波成分的幅值譜，而對(duì)于非 K 次諧波成分的頻譜則無(wú)法檢測(cè)出來(lái)，并且由于柵欄效應(yīng)和頻譜泄漏，對(duì)其它 K次諧波的幅值也有影響。 如將上述程序中的 N=100 改為 N=50，結(jié)果如圖 所示。 14 圖 N=50離散信號(hào)頻譜 可以看出，由于基頻頻率為 100/50=2Hz， x(t)中的 15Hz 分量未檢測(cè)出來(lái)，其它頻率分量的幅值也出現(xiàn)了誤差，這種柵欄效應(yīng)也是 DFT 應(yīng)用中不可避免的問(wèn)題之一。要減小柵欄效應(yīng)和頻譜泄漏，可以減小采樣頻率，增加采樣點(diǎn)數(shù)，以減小基頻值，使譜線變密，這樣原來(lái)漏掉的某些頻譜成分就可以檢測(cè)出來(lái)，但注意采樣頻率必須滿(mǎn)足采樣定理的要求，增加采樣點(diǎn)數(shù)也會(huì)增加系統(tǒng)計(jì)算 DFT 的開(kāi)銷(xiāo)。 本章小結(jié) 在本章中， 前五節(jié) 介紹了頻譜分析所用到的數(shù)學(xué)理論，分別為 采樣定理，快速傅立葉變換，離散傅立葉變換， 介紹了頻譜分析 基本原理以及 對(duì)數(shù)字信號(hào)處理的重要性， 闡述了頻譜分析的方法， 在最后的 種 對(duì)幾種 典型 信號(hào)進(jìn)行了簡(jiǎn)單的頻譜分析。 15 第 3 章 典型 信號(hào)去噪方法 IIR 和 FIR 濾波器的設(shè)計(jì) IIR 濾波器和 FIR 濾波器的設(shè)計(jì)方法完全不同。 IIR 濾波器設(shè)計(jì)方法有間接法和直接法，間接法是借助于模擬濾波器的設(shè)計(jì)方法進(jìn)行的。其設(shè)計(jì)步驟是：先設(shè)計(jì)過(guò)渡模擬濾波器得到系統(tǒng)函數(shù) H（ s），然后將 H（ s）按某種方法轉(zhuǎn)換成數(shù)字濾波器的系統(tǒng)函數(shù) H(z)。 FIR 濾波器比鞥采用間接法，常用的方法有窗函數(shù)法、頻率采樣發(fā)和切比雪夫等波紋逼近法 [14] 。對(duì)于線性相位濾波器，經(jīng)常采用 FIR 濾波器。 對(duì)于數(shù)字高通、帶通濾波器的設(shè)計(jì)，通用方法為雙線性變換法?？梢越柚谀M 濾 波器的頻率轉(zhuǎn) 換設(shè)計(jì)一個(gè)所需類(lèi)型的過(guò)渡模擬濾波器，再經(jīng)過(guò)雙線性變換將其轉(zhuǎn)換策劃那個(gè)所需的數(shù)字濾波器。具體設(shè)計(jì)步驟如下： （ 1）確定所需類(lèi)型數(shù)字濾波器的技術(shù)指標(biāo)。 （ 2）將所需類(lèi)型數(shù)字濾波器的邊界頻率轉(zhuǎn)換成相應(yīng)的模擬濾波器的邊界頻率，轉(zhuǎn)換公式為Ω =2/T tan( ) （ 3）將相應(yīng)類(lèi)型的模擬濾波器技術(shù)指標(biāo)轉(zhuǎn)換成模擬低通濾波器技術(shù)指標(biāo)。 （ 4）設(shè)計(jì)模擬低通濾波器。 （ 5）通過(guò)頻率變換將模擬低通轉(zhuǎn)換成相應(yīng)類(lèi)型的過(guò)渡模擬濾波器。 （ 6）采用雙線性變換法將相應(yīng)類(lèi)型的過(guò)渡模擬濾波器轉(zhuǎn)換成所需類(lèi)型的數(shù)字濾波器。 我們知道 ，脈沖響應(yīng)不變法的主要缺點(diǎn)是會(huì)產(chǎn)生頻譜混疊現(xiàn)象，使數(shù)字濾波器的頻響偏離模擬濾波器的頻響特性。為了克服之一缺點(diǎn)，可以采用雙線性變換法。 下面我們總結(jié)一下利用模擬濾波器設(shè)計(jì) IIR數(shù)字低通濾波器的步驟 [15] ： （ 1）確定數(shù)字低通濾波器的技術(shù)指標(biāo)：通帶邊界頻率、通帶最大衰減，阻帶截止頻 率、阻帶最小衰減。 （ 2）將數(shù)字低通濾波器的技術(shù)指標(biāo)轉(zhuǎn)換成相應(yīng)的模擬低通濾波器的技術(shù)指標(biāo)。 16 （ 3）按照模擬低通濾波器的技術(shù)指標(biāo)設(shè)計(jì)及過(guò)渡模擬低通濾波器。 （ 4）用雙線性變換法， 模擬濾波器系統(tǒng)函數(shù)轉(zhuǎn)換成數(shù)字低通濾波器系統(tǒng)函數(shù)。 如前所述， IIR 濾波器和 FIR 濾波器的設(shè)計(jì)方法有很大的區(qū)別。下面我們著重介紹用窗函數(shù)法設(shè)計(jì) FIR 濾波器的步驟。如下： （ 1）根據(jù)對(duì)阻帶衰減及過(guò)渡帶的指標(biāo)要求，選擇串窗數(shù)類(lèi)型（矩形窗、三角 窗 、漢寧窗、哈明窗、凱塞窗等），并估計(jì)窗口長(zhǎng)度 N。先按照阻帶衰減選擇窗函數(shù)類(lèi)型。原則是在保證阻帶衰減滿(mǎn)足要求的情況下，盡量選擇主瓣的窗函數(shù)。 （ 2）構(gòu)造希望逼近的頻率響應(yīng)函數(shù)。 （ 3）計(jì)算 h(n)。 （ 4）加窗得到設(shè)計(jì)結(jié)果。 在 Matlab 中，可以利用函數(shù) fir1 設(shè) 計(jì) FIR 濾波器，利用函數(shù) butter,cheby1和 ellip 設(shè)計(jì) IIR 濾波器，利用 Matlab 中的函數(shù) freqz 畫(huà)出各步步器的頻率響應(yīng)。 hn=fir1(M， wc， window)，可以指定窗函數(shù)向量 window。如果缺省 window參數(shù)，則 FIR 默認(rèn)為哈明窗。其中可選 的 Rectangular Barlrtt Hamming Hann Blackman 窗，其相應(yīng)的都有實(shí)現(xiàn)函數(shù)。 MATLAB 信號(hào)處理工具箱函數(shù) buttp buttor butter 是巴特沃斯濾波器設(shè)計(jì)函數(shù)，其有 5 種調(diào)用格式， 本文中 用到的 是[N,wc]=butter(N,wc,Rp,As,’ s’ ),該格式用于計(jì)算巴特沃斯模擬濾波器的階數(shù) N和 3dB 截止頻率 wc。 MATLAB 信號(hào)處理工具箱函數(shù) cheblap,cheblord 和 cheeby1 是切比雪夫 I 型濾波器設(shè)計(jì)函數(shù)。我們用到的是 cheeby1 函數(shù)，其調(diào)用格式如下：[B,A]=cheby1(N,Rp,wpo,’ ftypr’ ) [B,A]=cheby1(N,Rp,wpo,’ ftypr’ ,’ s’ ) 函數(shù) butter,cheby1 和 ellip 設(shè)計(jì) IIR濾波器時(shí)都是默認(rèn)的雙線性變換法，所以在設(shè)計(jì)濾波器時(shí) 只需要代入相應(yīng)的實(shí)現(xiàn)函數(shù)即可。 17 典型 信號(hào)的去噪對(duì)比 IIR 低通濾波器去噪 對(duì)正弦信號(hào)加噪聲之后 再 進(jìn)行去噪，利用 IIR濾波的方法 去除噪聲。 正弦信號(hào)為 x=A*cos(2*pi*f*t/fs)， 其中： A=1， f=1000， fs=10^5。 程序代碼如下： clf。 Ft=9000。 Fp=100。 Fs=300。 wp=2*pi*Fp/Ft。 ws=2*pi*Fs/Ft。 fp=2*Ft*tan(wp/2)。 fs=2*Fs*tan(wp/2)。 [n11,wn11]=buttord(wp,ws,1,50,39。s39。)。 %求低通濾波器的階數(shù)和截止頻率 [b11,a11]=butter(n11,wn11,39。s39。)。 %求 S域的頻率響應(yīng)的參數(shù) [num11,den11]=bilinear(b11,a11,)。 %雙線性變換實(shí)現(xiàn) S域到 Z域的變換 [h,w]=freqz(num11,den11)。 %根據(jù)參數(shù)求出頻率響應(yīng) plot(w*8000*,abs(h))。 legend(39。用 butter設(shè)計(jì) 39。)。 grid。 A=1。%信號(hào)的幅值 f=1000。%信號(hào)的頻率 fs=10^5。 %采樣頻率 t=(0:999)。 %采樣點(diǎn) 18 Mlag=100。 %相關(guān)函數(shù)長(zhǎng)度變量 x=A*cos(2*pi*f*t/fs)。 %輸入正弦波信號(hào) xmean=mean(x)。 %正弦波信號(hào)均值 xvar=var(x,1)。%正弦波信號(hào)方差 xn=awgn(x,5)。 %給正弦波信號(hào)加入信噪比為 20dB的高斯白噪聲 figure(1) plot(t,xn)%繪制輸入信號(hào)圖像 title(39。輸入信號(hào)圖像 39。) xlabel(39。x軸單位： t/s39。,39。color39。,39。b39。) ylabel(39。y軸單位： f/HZ39。,39。color39。,39。b39。) S=fft(x)。 %傅里葉變換 z11=filter(num11,den11,xn)。 m11=fft(z11)。 %求濾波后的信號(hào) figure(1)。 subplot(2,1,1)。 plot(abs(S),39。g39。)。 title(39。濾波前信號(hào)的頻譜 39。)。 grid。 figure(1)。 subplot(2,1,2)。 plot(abs(m11),39。r39。)。 title(39。濾波后信號(hào)的頻譜 39。)。 grid。 figure(2)。 19 subplot(2,1,1)。 plot(xn)。 title(39。濾波前信號(hào)的波形 39。)。 grid。 figure(2)。 subplot(2,1,2)。 plot(z11)。 title(39。濾波后的信號(hào)波形 39。)。 0 100 200 300 400 500 600 700 800 900 10000200400600濾波前信號(hào)的頻譜0 100 200 300 400 500 600 700 800 900 10000200400600濾波后信號(hào)的頻譜 圖 正弦 信號(hào)頻譜圖 圖 為濾波 前 后的頻譜對(duì)比圖。 20 0 100 200 300 400 500 600 700 800 900 100042024濾波前信號(hào)的波形0 100 200 300 400 500 600 700 800 900 100021012濾波后的信號(hào)波形 圖 去噪前后噪聲對(duì)比 圖 圖 為濾波前后的信號(hào) 波形 對(duì)比圖，可以看到，利用 IIR 低通濾波器雖然仍有小部分的噪聲，信號(hào)并不完全平滑，但是絕大部分的 噪聲已經(jīng)被濾除。 FIR 低通濾波器去噪 設(shè)計(jì) FIR 濾波器對(duì)信號(hào)進(jìn)行去噪，信號(hào)為加噪聲的正弦信號(hào)，信號(hào)表達(dá)式為x=cos(2*pi*f*t/fs)，其中， f=1000， fs=10^5。 程序代碼如下： A=1。%信號(hào)的幅值 f=1000。%信號(hào)的頻率 fs=10^5。 %采樣頻率 t=(0:999)。 %采樣點(diǎn) Mlag=100。 %相關(guān)函數(shù)長(zhǎng)度變量 x=A*cos(2*pi*f*t/fs)。 %輸入正弦波信號(hào) xmean=mean(x)。 %正弦波信號(hào)均值 21 xvar=var(x,1)。%正弦波信號(hào)方差 xn=awgn(x,5)。 %給正弦波信號(hào)加入信噪比為 20dB的高斯白噪聲 figure(1) plot(t,xn)%繪制輸入信號(hào)圖像 title(39。輸入信號(hào)圖像 39。) xlabel(39。x軸單位： t/s39。,39。color39。,39。b39。) ylabel(39。y軸單位： f/HZ39。,39。color39。,39。b39。) %fir濾波 wp=*pi。 %通帶截止頻率 ws=*pi。 %阻帶截止頻率 DB=wswp。 %過(guò)渡帶寬度 N0=ceil(*pi/DB)。 M=N0+mod(N0+1,2)。%計(jì)算 fir濾波器階數(shù) wc=(wp+ws)/2/pi。%計(jì)算理想低通濾波器通帶截止頻率（關(guān)于π歸一化） hn=fir1(M,wc)。%調(diào)用 fir1計(jì)算 FIRDF的 h(n) y1n=filter(hn,1,xn)。 %將輸入信號(hào)通過(guò) fir濾波器 figur