喬建華+張雪英
摘 要: 譜分析是數字信號處理中的一個重要問題,初學者普遍對連續信號譜分析理解不深,尤其是在誤差分析時缺乏統一示例,更容易產生困惑。介紹了用離散傅里葉變換(DFT)對連續信號進行譜分析的過程,并詳細說明了誤差產生的原因和減小誤差的方法。而且通過對模擬信號譜分析的實例全面說明了各項誤差的影響及解決方案,并應用Matlab直觀地進行了分析和對比驗證。
關鍵詞: 譜分析; 誤差分析; DFT; Matlab
中圖分類號: TN911.7?34 文獻標識碼: A 文章編號: 1004?373X(2014)13?0053?04
Error of continuous signal spectrum analysis based on DFT
QIAO Jian?hua1,2, ZHANG Xue?ying2
(1. College of Electronic & Information Engineering, Taiyuan University of Science & Technology, Taiyuan 030024, China;
2. College of Information Engineering, Taiyuan University of Technology, Taiyuan 030024, China)
Abstract: Spectrum analysis is one of the important missions in digital signal processing. The understanding of most beginners is not deep to continuous signal spectrum analysis, especially lacking of a uniform example in the error analysis, so it is easier to be confused. The process of continuous time signal spectrum analysis using the discrete Fourier transform (DFT) is introduced in this paper. The reasons of generating the errors and the methods of minimizing the errors are illustrated in detail. The influence and solutions of various errors are elaborated through a concrete example of analog signal spectrum analysis. It is verified intuitively through analysis and contrast with Matlab software.
Keywords: spectrum analysis; analysis error; DFT; Matlab
0 引 言
連續信號的譜分析是數字信號處理的一個重要內容。對信號進行譜分析,關鍵是得到信號的傅里葉變換。由于離散傅里葉變換(DFT)有快速算法FFT,因此經常用DFT(FFT)對信號進行譜分析。但是根據傅里葉變換理論,若信號持續時間有限長,則其頻譜無限寬;若信號的頻譜有限寬,則其持續時間無限長[1]。因此,用有限長序列的有限點的DFT對連續信號進行譜分析必然是近似的,要明確誤差產生的原因,并能夠采取合理措施來降低誤差。文獻[2?5]也以不同的側重點介紹了DFT在譜分析中的仿真實現,物理實現,具體應用以及針對各種不同信號的譜分析。而對于實際中廣泛存在的連續信號,初學者普遍對其應用DFT的譜分析的原理和過程理解不深,尤其對誤差產生原因不明確,對減小誤差方案不清晰,在具體的應用中不知如何整合這些方案,這都是對整個譜分析過程研究不透徹的緣故。本文從具體的模擬信號譜分析的實例入手,來充分說明這一點,并通過Matlab演示來加深對整個譜分析過程中誤差的認識。
1 用DFT對連續信號譜分析
假設模擬信號為[xa(t)],根據Nyquist采樣定理,使用模擬信號帶寬2倍以上的采樣率才能保證頻譜不混疊,但實際上理想的帶限信號時域上總是無限長的,因此在工程應用時對于帶寬外能量占總能量比例很小的信號,可以近似看成是帶限信號,認為抽樣信號是不失真的。
例如設[xa(t)=e-1 000|t|,]其傅里葉變換為:
[Xa(jΩ)=FT[xa(t)]=-∞∞xa(t)e-jΩtdt=-∞∞e-1 000|t|e-jΩtdt=2 0001 0002+Ω2] (1)
可以看到這個信號不是一個理想帶限的信號。但是隨著[Ω]增加,[Xa(jΩ)]會逐漸減小,可以規定一個帶寬[Ωc,]使得[Xa(jΩc)?1,]認為當采樣頻率[fs>Ωcπ]時可以無失真地恢復原模擬信號,如本例中取[fc=Ωc2π=2]kHz,此即造成了誤差的產生。[xa(t)]及其幅頻特性曲線如圖1(a),(b)所示。
對[xa(t)]以采樣頻率[fs=5 ]kHz采樣,即采樣間隔[T=][0.2 ]ms,得采樣信號[xa(nT)=xa(t)t=nT=x(n),]在式(1)中,由于[t→nT,dt→T,-∞∞dt→n=-∞∞T,]因此得采樣信號的頻譜:
[Xa(jΩ)≈T?n=-∞∞xa(nT)e-jΩnT] (2)
假設將采樣序列[xa(nT)]截取成長度為[M]的有限長序列,此截斷即為誤差產生之二。若頻率用[f]表示,則式(2)寫為:
[Xa(jf)≈T?n=0M-1xa(nT)e-j2πfnT] (3)
顯然[Xa(jf)]仍是[f]的連續周期函數,[xa(nT)]和[Xa(jf)]如圖1(c),(d)所示。
為了進行數值運算,頻域上也要離散化,因此在頻域的一個周期[0, fs]內等間隔采樣[N]點,如圖1(f)所示。由于間隔采樣有可能漏掉重要的頻域信息,成為誤差產生之三。設頻域采樣間隔為[F,]則[F=fsN=1(NT),]將[f=kF]代入式(3)中可得對采樣信號頻譜[Xa(j f)]的采樣:
[Xa(jkF)≈T?n=0M-1xa(nT)e-j2πNkn] (4)
根據頻域采樣定理[6],如果序列[x(n)]的長度為[M,]只有當頻域采樣點數[N≥M]時,才可由頻域采樣[X(k)]恢復原序列[x(n),]否則會產生時域混疊。因此為簡便計算,將序列也以[N]點長度來截取,則式(4) 的求和上限為[N-1,]并令:[Xa(k)=Xa(jkF),x(n)=xa(nT)],則:
[Xa(k)≈T?n=0N-1x(n)e-j2πNkn=T?DFT[x(n)]N, k=0,1,…,N-1] (5)
可見,對近似持續時間有限的帶限連續信號,其頻譜特性可以通過對連續信號采樣[N]點并進行DFT再乘以[T]的近似方法得到。因此用DFT對模擬信號譜分析,產生誤差是必然的。其具體分析過程是首先確定模擬信號的最高頻率[fc]和頻率分辨率[F,]再選取參數,包括采樣頻率[fs,]采樣點數[N,]記錄時間[Tp]等。其取值分別為:
[fs≥2fcN=fsF≥2fcFTp=NT=Nfs≥1F] (6)
然后根據計算的參數,通過采樣得到信號序列[x(n)],即可根據式(5)進行運算分析。
圖1 DFT對連續信號譜分析過程
下面具體分析誤差的影響及其減小誤差的措施。
2 誤差分析
2.1 頻譜混疊
在對模擬信號[xa(t)]進行采樣時,必須滿足采樣定理,即采樣頻率[fs≥2fc,]而通過上面分析時域有限的信號不可能是銳截止的,并且信號中不可避免地有一些高頻雜散信號,因此在采樣之前,一般都要對模擬信號進行濾波,濾除高頻雜散信號。如本例中,[fc]取2 kHz,當[fs]=5 kHz時,采樣信號的頻譜混疊很小,見圖1(d);但當選[fs]=1 kHz時,其頻譜如圖1(h)所示,混疊嚴重。以折疊頻率[fs2]處的頻譜幅度來比較[T=]0.2 ms和[T=]1 ms時的混疊程度,此處[k=256,]鍵入Matlab語句,得到:
abs(X1(k))/max(abs(X1))
ans=0.009 4
abs(X2(k))/max(abs(X2))
ans=0.212 6
而在[fc]處,模擬信號的幅度比值為:
abs(Xa(2k))/max(abs(Xa))
ans=0.006 6
可以看出,隨著采樣頻率的減小,混疊現象加大。
因此,要減小混疊,必須滿足乃奎斯特采樣定理,并且在采樣前進行預濾波,濾除高于折疊頻率[fs2]的頻率成分,一般取采樣頻率[6][fs≥(3~5)fc。]
2.2 截斷效應
對無限長的模擬信號,用DFT進行譜分析時,必須先進行截斷,通過采樣才能得到有限點的序列,這樣必然產生誤差。截斷可以理解為加窗,即:
[y(n)=x(n)w(n)] (7)
式中:[x(n)]為模擬信號經采樣得到的時域離散信號;[w(n)]為窗函數序列。根據頻域卷積定理,加窗后的信號頻譜為:
[Y(ejω)=12πX(ejω)*W(ejω)] (8)
顯然與原序列的頻譜是不同的。
文獻[7]對常用的五種窗函數的時頻域波形有詳細的對比,是對窗函數本身的性能做了說明。下面僅以矩形窗和hamming窗為例來從加窗對信號頻譜的影響方面來進行分析。圖2所示為矩形窗和hamming窗的時域和頻域波形。由時域波形可見,矩形窗是對原信號的原樣截取,hamming窗是在截取的同時對信號做了加權。由頻域波形可見,矩形窗的主瓣窄,其寬度[8]為[4πN,]但旁瓣幅度也高;hamming窗是以主瓣展寬為代價換取了旁瓣幅度的減小,其主瓣寬度[8]為[8πN,]而主瓣寬度影響的就是頻率分辨率。可見,通過加窗對信號截斷使得對模擬信號的譜分析產生誤差,表現為原本離散的譜線展寬(高頻泄漏)和譜間干擾等現象,即所謂截斷效應。從而降低了譜的分辨率,使頻率相近的兩個信號不易分清。
圖2 矩形窗和hamming窗的時頻曲線
因此,為了減小截斷效應的影響,一種方法是增大截取長度,通過改變[N,]使窗函數的主瓣變窄,提高頻率分辨率。另一種方法是改變截斷窗函數的形狀。
下面通過一個例子,來更清楚地展示泄漏和譜間干擾的影響,以及通過截取長度和加不同的窗函數對頻譜的改善作用[9],設:
[x(n)=cos(2πf1n)+sin(2πf2n)+cos(2πf3n)]
其中[f1]=25 Hz,[f2]=50 Hz,[f3]=100 Hz,可見其最高頻率為[f3,]頻率分辨率至少為25 Hz,則根據式(6),截取長度[Tp≥1F=]0.04 s,取采樣頻率[fs]=400 Hz,則最小采樣點[N=][Tp?][fs]=16點,圖3所示為其加矩形窗和hamming窗以[N,][4N,][8N]點分別截斷的頻譜效果。由圖可以清晰地看到,當截取長度太短的時候,泄漏和譜間干擾都非常嚴重,25 Hz和50 Hz的兩條譜線已無法分辨,只有增加截取長度[N,]才使得主瓣變窄,提高了頻率分辨率;而采用hamming窗又降低了旁瓣的影響,減小了譜間干擾。而且可以看出,當[N]一定時,旁瓣越小的窗函數,其主瓣就越寬。增加[N]使主瓣變窄,但旁瓣的相對幅度并不減小。矩形窗的主瓣較窄,但是旁瓣幅度也較高。可見,當截取長度一定時,頻率分辨率和譜間干擾是相互抵消的,只能以一種優勢來換取另一種性能的降低。
圖3 矩形窗和hamming窗以不同長度截斷效果
2.3 柵欄效應
柵欄效應是指用DFT對連續信號進行譜分析時,由于DFT的離散特性,使離散點之間的頻譜無法得到,相當于透過柵欄觀察頻譜,只看到等間隔的離散點的頻線,其他的頻譜都被柵欄擋住了,故稱之為柵欄效應。因此用DFT得到的離散譜線的包絡只能是近似譜。
為了減小柵欄效應,可以多增加些柵欄的縫隙,即增大DFT變換點數,這一方面可以通過在原序列尾部補零來實現[10]。例如對下式信號進行分析:
[x(n)=sin(2πf1n)+sin(2πf2n)+sin(2πf3n)]
其中[f1]=20 Hz,[f2]=20.5 Hz,[f3]=40 Hz,由其最高頻率為[f3],頻率分辨率至少要達到0.5 Hz,根據式(6)選擇參數,若采樣頻率[fs]取80 Hz,則采樣點數[N≥][800.5=]160。若[N]取160點,再進行DFT,可得到圖4波形[x1(n)]及幅度譜[X1(k),]由[X1(k)]曲線可見,無法分辨20 Hz和20.5 Hz的兩條譜線,這就是由于柵欄效應造成的。如果在序列后補零加長到512點,如圖4中[x2(n),]再看其幅度譜[X2(k),]頻譜被細化,挨近的兩條譜線已能有所區分。
圖4 補零抑制柵欄效應
但是,如果對原序列采樣128點,即使補零加長到512點,如圖5所示,其頻譜[X2(k)]也無法分辨出這兩條譜線。這是因為截斷效應已經使得頻譜變模糊,即使再增加柵欄的縫隙,看到的也是模糊的頻譜,補零并沒有真正提高頻率分辨率。這種情況下,只能通過增加截取長度,首先提高頻譜分辨率。如對信號直接采樣512點,觀察圖5的[X3(k)]幅頻特性曲線,此時,這兩條譜線可以清晰地分辨出來。
圖5 增加截取長度抑制柵欄效應
總之,對模擬信號的DFT譜分析主要集中于兩個方面:首先是對信號的采樣,其中涉及的問題有頻譜混疊和截斷效應,因此要選擇合適的采樣頻率和窗函數,足夠的記錄時間;然后是對得到的時域離散信號[x(n)]進行DFT運算,這時要注意柵欄效應,要計算足夠點數的DFT,避免漏掉鄰近的譜線。因此在選擇參數的時候,在依據式(6)的基礎上,實際應用中要適當加大取值。
3 結 語
本文針對數字信號處理中的一個普遍應用——采用DFT對連續信號譜分析中的誤差問題,從基本原理上說明了利用DFT對連續信號譜分析的近似性,對造成誤差的三種主要因素(頻譜混疊、截斷效應和柵欄效應)的產生做了深入分析,提出減小誤差的相應措施,以及參數選擇的基本原則,而且以詳實細致的圖例對相應問題作了說明。需要提出的是如果利用近代譜估計的方法,可以得到一些更好的效果,下一階段可以在這方面繼續進行研究。
參考文獻
[1] OPPENHEIM A V, SCHAFER R W.數字信號處理[M].董士嘉,楊耀增,譯.北京:科學出版社,1980.
[2] 謝海霞,孫志雄.連續非周期信號頻譜分析及Matlab實現[J].現代電子技術,2013,36(11):53?56.
[3] 郭連平,田書林,王志剛.數字信號處理過程中信號截位誤差抑制方法研究[J].信號處理,2013(5):544?546.
[4] 溫萬惠.在生理信號頻譜分析中掌握離散傅里葉變換法[J].西南師范大學學報:自然科學版,2010(6):227?230.
[5] 林愛英,滕紅麗,袁超,等.DFT在信號譜分析中的應用[J].安徽工業大學學報:自然科學版,2011(2):192?196.
[6] 高西全,丁玉美.數字信號處理[M].3版.西安:西安電子科技大學出版社,2008.
[7] 汪偉,謝皓臣,梁光明,等.加窗離散傅里葉變換性能分析和比對[J].現代電子技術,2012,35(3):115?118.
[8] ORFANIDIS S J.信號處理導論(影印版)[M].北京:清華大學出版社,1999.
[9] 陳懷琛.Matlab及在電子信息課程中的應用[M].3版.北京:電子工業出版社,2006.
[10] 薛年喜.Matlab在數字信號處理中的應用[M].2版.北京:清華大學出版社,2008.
圖3 矩形窗和hamming窗以不同長度截斷效果
2.3 柵欄效應
柵欄效應是指用DFT對連續信號進行譜分析時,由于DFT的離散特性,使離散點之間的頻譜無法得到,相當于透過柵欄觀察頻譜,只看到等間隔的離散點的頻線,其他的頻譜都被柵欄擋住了,故稱之為柵欄效應。因此用DFT得到的離散譜線的包絡只能是近似譜。
為了減小柵欄效應,可以多增加些柵欄的縫隙,即增大DFT變換點數,這一方面可以通過在原序列尾部補零來實現[10]。例如對下式信號進行分析:
[x(n)=sin(2πf1n)+sin(2πf2n)+sin(2πf3n)]
其中[f1]=20 Hz,[f2]=20.5 Hz,[f3]=40 Hz,由其最高頻率為[f3],頻率分辨率至少要達到0.5 Hz,根據式(6)選擇參數,若采樣頻率[fs]取80 Hz,則采樣點數[N≥][800.5=]160。若[N]取160點,再進行DFT,可得到圖4波形[x1(n)]及幅度譜[X1(k),]由[X1(k)]曲線可見,無法分辨20 Hz和20.5 Hz的兩條譜線,這就是由于柵欄效應造成的。如果在序列后補零加長到512點,如圖4中[x2(n),]再看其幅度譜[X2(k),]頻譜被細化,挨近的兩條譜線已能有所區分。
圖4 補零抑制柵欄效應
但是,如果對原序列采樣128點,即使補零加長到512點,如圖5所示,其頻譜[X2(k)]也無法分辨出這兩條譜線。這是因為截斷效應已經使得頻譜變模糊,即使再增加柵欄的縫隙,看到的也是模糊的頻譜,補零并沒有真正提高頻率分辨率。這種情況下,只能通過增加截取長度,首先提高頻譜分辨率。如對信號直接采樣512點,觀察圖5的[X3(k)]幅頻特性曲線,此時,這兩條譜線可以清晰地分辨出來。
圖5 增加截取長度抑制柵欄效應
總之,對模擬信號的DFT譜分析主要集中于兩個方面:首先是對信號的采樣,其中涉及的問題有頻譜混疊和截斷效應,因此要選擇合適的采樣頻率和窗函數,足夠的記錄時間;然后是對得到的時域離散信號[x(n)]進行DFT運算,這時要注意柵欄效應,要計算足夠點數的DFT,避免漏掉鄰近的譜線。因此在選擇參數的時候,在依據式(6)的基礎上,實際應用中要適當加大取值。
3 結 語
本文針對數字信號處理中的一個普遍應用——采用DFT對連續信號譜分析中的誤差問題,從基本原理上說明了利用DFT對連續信號譜分析的近似性,對造成誤差的三種主要因素(頻譜混疊、截斷效應和柵欄效應)的產生做了深入分析,提出減小誤差的相應措施,以及參數選擇的基本原則,而且以詳實細致的圖例對相應問題作了說明。需要提出的是如果利用近代譜估計的方法,可以得到一些更好的效果,下一階段可以在這方面繼續進行研究。
參考文獻
[1] OPPENHEIM A V, SCHAFER R W.數字信號處理[M].董士嘉,楊耀增,譯.北京:科學出版社,1980.
[2] 謝海霞,孫志雄.連續非周期信號頻譜分析及Matlab實現[J].現代電子技術,2013,36(11):53?56.
[3] 郭連平,田書林,王志剛.數字信號處理過程中信號截位誤差抑制方法研究[J].信號處理,2013(5):544?546.
[4] 溫萬惠.在生理信號頻譜分析中掌握離散傅里葉變換法[J].西南師范大學學報:自然科學版,2010(6):227?230.
[5] 林愛英,滕紅麗,袁超,等.DFT在信號譜分析中的應用[J].安徽工業大學學報:自然科學版,2011(2):192?196.
[6] 高西全,丁玉美.數字信號處理[M].3版.西安:西安電子科技大學出版社,2008.
[7] 汪偉,謝皓臣,梁光明,等.加窗離散傅里葉變換性能分析和比對[J].現代電子技術,2012,35(3):115?118.
[8] ORFANIDIS S J.信號處理導論(影印版)[M].北京:清華大學出版社,1999.
[9] 陳懷琛.Matlab及在電子信息課程中的應用[M].3版.北京:電子工業出版社,2006.
[10] 薛年喜.Matlab在數字信號處理中的應用[M].2版.北京:清華大學出版社,2008.
圖3 矩形窗和hamming窗以不同長度截斷效果
2.3 柵欄效應
柵欄效應是指用DFT對連續信號進行譜分析時,由于DFT的離散特性,使離散點之間的頻譜無法得到,相當于透過柵欄觀察頻譜,只看到等間隔的離散點的頻線,其他的頻譜都被柵欄擋住了,故稱之為柵欄效應。因此用DFT得到的離散譜線的包絡只能是近似譜。
為了減小柵欄效應,可以多增加些柵欄的縫隙,即增大DFT變換點數,這一方面可以通過在原序列尾部補零來實現[10]。例如對下式信號進行分析:
[x(n)=sin(2πf1n)+sin(2πf2n)+sin(2πf3n)]
其中[f1]=20 Hz,[f2]=20.5 Hz,[f3]=40 Hz,由其最高頻率為[f3],頻率分辨率至少要達到0.5 Hz,根據式(6)選擇參數,若采樣頻率[fs]取80 Hz,則采樣點數[N≥][800.5=]160。若[N]取160點,再進行DFT,可得到圖4波形[x1(n)]及幅度譜[X1(k),]由[X1(k)]曲線可見,無法分辨20 Hz和20.5 Hz的兩條譜線,這就是由于柵欄效應造成的。如果在序列后補零加長到512點,如圖4中[x2(n),]再看其幅度譜[X2(k),]頻譜被細化,挨近的兩條譜線已能有所區分。
圖4 補零抑制柵欄效應
但是,如果對原序列采樣128點,即使補零加長到512點,如圖5所示,其頻譜[X2(k)]也無法分辨出這兩條譜線。這是因為截斷效應已經使得頻譜變模糊,即使再增加柵欄的縫隙,看到的也是模糊的頻譜,補零并沒有真正提高頻率分辨率。這種情況下,只能通過增加截取長度,首先提高頻譜分辨率。如對信號直接采樣512點,觀察圖5的[X3(k)]幅頻特性曲線,此時,這兩條譜線可以清晰地分辨出來。
圖5 增加截取長度抑制柵欄效應
總之,對模擬信號的DFT譜分析主要集中于兩個方面:首先是對信號的采樣,其中涉及的問題有頻譜混疊和截斷效應,因此要選擇合適的采樣頻率和窗函數,足夠的記錄時間;然后是對得到的時域離散信號[x(n)]進行DFT運算,這時要注意柵欄效應,要計算足夠點數的DFT,避免漏掉鄰近的譜線。因此在選擇參數的時候,在依據式(6)的基礎上,實際應用中要適當加大取值。
3 結 語
本文針對數字信號處理中的一個普遍應用——采用DFT對連續信號譜分析中的誤差問題,從基本原理上說明了利用DFT對連續信號譜分析的近似性,對造成誤差的三種主要因素(頻譜混疊、截斷效應和柵欄效應)的產生做了深入分析,提出減小誤差的相應措施,以及參數選擇的基本原則,而且以詳實細致的圖例對相應問題作了說明。需要提出的是如果利用近代譜估計的方法,可以得到一些更好的效果,下一階段可以在這方面繼續進行研究。
參考文獻
[1] OPPENHEIM A V, SCHAFER R W.數字信號處理[M].董士嘉,楊耀增,譯.北京:科學出版社,1980.
[2] 謝海霞,孫志雄.連續非周期信號頻譜分析及Matlab實現[J].現代電子技術,2013,36(11):53?56.
[3] 郭連平,田書林,王志剛.數字信號處理過程中信號截位誤差抑制方法研究[J].信號處理,2013(5):544?546.
[4] 溫萬惠.在生理信號頻譜分析中掌握離散傅里葉變換法[J].西南師范大學學報:自然科學版,2010(6):227?230.
[5] 林愛英,滕紅麗,袁超,等.DFT在信號譜分析中的應用[J].安徽工業大學學報:自然科學版,2011(2):192?196.
[6] 高西全,丁玉美.數字信號處理[M].3版.西安:西安電子科技大學出版社,2008.
[7] 汪偉,謝皓臣,梁光明,等.加窗離散傅里葉變換性能分析和比對[J].現代電子技術,2012,35(3):115?118.
[8] ORFANIDIS S J.信號處理導論(影印版)[M].北京:清華大學出版社,1999.
[9] 陳懷琛.Matlab及在電子信息課程中的應用[M].3版.北京:電子工業出版社,2006.
[10] 薛年喜.Matlab在數字信號處理中的應用[M].2版.北京:清華大學出版社,2008.