一、DTFT、DFT、FFT、DFS的關系
1、DTFT與DFT的關系
在實際應用中,計算機隻能處理離散信号,是以對連續信号 x ( t ) x(t) x(t)進行時域采樣,得到一組離散樣本 x ( n ) x(n) x(n),對它進行傅裡葉變換得到
X ( w ) = ∑ n = − ∞ ∞ x ( n ) e − j w n X(w)=\sum_{n=-∞}^{∞} x(n)e^{-jwn} X(w)=n=−∞∑∞x(n)e−jwn
上式即為離散時間傅裡葉變換(DTFT),由于變換後得到的頻域值仍然是連續的,繼續對頻域進行采樣,得到
X ( k ) = ∑ n = 0 N − 1 x ( n ) e − j 2 π k n N X(k)=\sum_{n=0}^{N-1} x(n)e^{-j\frac{2\pi k n}{N}} X(k)=n=0∑N−1x(n)e−jN2πkn
上式就是離散傅裡葉變換(DFT)。目前計算機中常用的快速傅裡葉變換(FFT),就是DFT的快速算法。
也就是說:DFT是DTFT的離散化,DTFT在頻域是連續的,在歸一化ω軸上是連續的,以2π為周期。
2、DFS與DFT的關系
離散傅裡葉級數(DFS)是針對周期序列的,即無限長序列;
而将DFS取出一個周期就是DFT,即DFT是有限長序列;
實際上,DFT就是DFS在一個周期内的取值,一個序列的DFT實際上就是這個序列的頻譜在一個周期内等間隔采樣的樣點值。
二、DFT的實質
離散傅裡葉變換:DFT(Discrete Fourier Transform)
其實質是有限長序列傅裡葉變換的有限點離散采樣,進而實作了頻域離散化,使數字信号處理可以在頻域采用數值運算的方法進行,這樣就大大增加了數字信号處理的靈活性。
更重要的是,DFT有多種快速算法,統稱為快速傅裡葉變換(Fast Fourier Transform,FFT),進而使信号的實時處理和裝置的簡化得以實作。
三、DFT的定義
設 x ( n ) x(n) x(n)是一個長度為 M M M的有限長序列,則定義 x ( n ) x(n) x(n)的 N N N點離散傅裡葉變換為
X ( k ) = D F T [ x ( n ) ] = ∑ n = 0 N − 1 x ( n ) W N k n k = 0 , 1 , ⋯ , N − 1 X(k)=DFT[x(n)]=\sum_{n=0}^{N-1}x(n)W_{N}^{kn} \qquad k=0,1,\cdots,N-1 X(k)=DFT[x(n)]=n=0∑N−1x(n)WNknk=0,1,⋯,N−1
X ( k ) X(k) X(k)的離散傅裡葉逆變換(Inverse Discrete Fourier Transfrom,IDFT)為
x ( n ) = I D F T [ X ( k ) ] = 1 N ∑ k = 0 N − 1 X ( k ) W N − k n n = 0 , 1 , ⋯ , N − 1 x(n)=IDFT[X(k)]={1\over N}\sum_{k=0}^{N-1}X(k)W_{N}^{-kn} \qquad n=0,1,\cdots,N-1 x(n)=IDFT[X(k)]=N1k=0∑N−1X(k)WN−knn=0,1,⋯,N−1
式中, W N = e − j 2 π N W_{N}=e^{-j \frac{2\pi}{N}} WN=e−jN2π, N N N稱為 D F T DFT DFT變換區間長度, N ≥ M N\ge M N≥M。
四、DFT的實體意義
DFT與傅裡葉變換和Z變換的關系:
設序列 x ( n ) x(n) x(n)的長度為 M M M,其 Z Z Z變換和 N ( N ≥ M ) N(N\ge M) N(N≥M)點DFT分别為
X ( z ) = Z T [ x ( n ) ] = ∑ n = 0 M − 1 x ( n ) z − n X(z)=ZT[x(n)]=\sum_{n=0}^{M-1}x(n)z^{-n} X(z)=ZT[x(n)]=n=0∑M−1x(n)z−n
X ( k ) = D F T [ x ( n ) ] N = ∑ n = 0 M − 1 x ( n ) W N k n k = 0 , 1 , ⋯ , N − 1 X(k)=DFT[x(n)]_N=\sum_{n=0}^{M-1}x(n)W_{N}^{kn} \qquad k=0,1,\cdots,N-1 X(k)=DFT[x(n)]N=n=0∑M−1x(n)WNknk=0,1,⋯,N−1
比較上面二式可得關系式
X ( k ) = X ( z ) ∣ z = e j 2 π N k k = 0 , 1 , ⋯ , N − 1 (1) X(k)=X(z)\Big|_{z=e^{j\frac{2\pi}{N}k}} \qquad k=0,1,\cdots,N-1 \tag{1} X(k)=X(z)∣∣∣z=ejN2πkk=0,1,⋯,N−1(1)
或
X ( k ) = X ( e j ω ) ∣ ω = 2 π N k k = 0 , 1 , ⋯ , N − 1 (2) X(k)=X(e^{j\omega})\Big|_{\omega=\frac{2\pi}{N}k} \qquad k=0,1,\cdots,N-1 \tag{2} X(k)=X(ejω)∣∣∣ω=N2πkk=0,1,⋯,N−1(2)
(1)式表明序列 x ( n ) x(n) x(n)的 N N N點 D F T DFT DFT是 x ( n ) x(n) x(n)的 Z Z Z變換在機關圓上的 N N N點等間隔采樣。(2)式則說明 X ( k ) X(k) X(k)為 x ( n ) x(n) x(n)的傅裡葉變換 X ( e j ω ) X(e^{j\omega}) X(ejω)在區間 [ 0 , 2 π ] [0,2\pi] [0,2π]上的 N N N點等間隔采樣。這就是DFT的實體意義。
由此可見,DFT的變化區間長度N不同,表示對 X ( e j ω ) X(e^{j\omega}) X(ejω)在區間 [ 0 , 2 π ] [0,2\pi] [0,2π]上的采樣間隔和采樣點數不同,是以DFT的變換結果不同。
xn=boxcar(4);
M=1024; % DFT變換區間長度
Xw=fft(xn,M); % fft函數在xn後面補零,傳回xn的M點DFT變換結果向量
wk=2*(0:M-1)/M;
X8k=fft(xn,8); % 計算xn的8點DFT
X16k=fft(xn,16); % 計算xn的16點DFT
subplot(3,2,1);
plot(wk,abs(Xw));
title('(a)x(n)的幅頻特性曲線 ');xlabel('\omega/\pi');ylabel('|X(e^j^\omega)|');axis([0,2,0,4.5]);
subplot(3,2,3);
k=0:7; stem(k,abs(X8k),'.');
title('(b)x(n)的8點DFT');xlabel('k');ylabel('|X(k)|');axis([0,8,0,4.5]);
subplot(3,2,5);
k=0:15; stem(k,abs(X16k),'.')
title('(c)x(n)的16點DFT');xlabel('k');ylabel('|X(k)|');axis([0,16,0,4.5]);
%%%%%%%%%% DFT的MATLAB計算示例 %%%%%%%%%%
xn=[1 1 1 1]; %輸入長度M=4的時域序列向量
Xk16=fft(xn,16); %計算xn的16點DFT
Xk32=fft(xn,32); %計算xn的32點DFT
k=0:15;
wk=2*k/16; %産生16點DFT對應的采樣點頻率(關于π歸一化值)
subplot(2,2,1);
stem(wk,abs(Xk16),'.'); %繪制16點DFT的幅頻特性圖
title('(a)16點DFT的幅頻特性圖');
xlabel('ω/π');
ylabel('幅度');
subplot(2,2,3);
stem(wk,angle(Xk16),'.'); %繪制16點DFT的相頻特性圖
line([0,2],[0,0]); %繪制一條從(0,0)到(2,0)的線條
title('(b)16點DFT的相頻特性圖')
xlabel('ω/π');
ylabel('相位');
axis([0,2,-3.5,3.5]); %更改坐标軸範圍,使 x 軸的範圍從0到2,y軸的範圍從-3.5 到3.5
k=0:31;
wk=2*k/32; %産生32點DFT對應的采樣點頻率(關于π歸一化值)
subplot(2,2,2);
stem(wk,abs(Xk32),'.'); %繪制32點DFT的幅頻特性圖
title('(c)32點DFT的幅頻特性圖');
xlabel('ω/π');
ylabel('幅度');
subplot(2,2,4);
stem(wk,angle(Xk32),'.'); %繪制32點DFT的相頻特性圖
line([0,2],[0,0]);
title('(d)32點DFT的相頻特性圖');
xlabel('ω/π');
ylabel('相位');
axis([0,2,-3.5,3.5]);
五、用DFT對信号進行譜分析
所謂信号的譜分析,就是計算信号的傅裡葉變換。連續信号與系統的傅裡葉分析不便于直接用計算機進行計算,使其應用受到限制。而DFT是一種時域和頻域均離散化的變換,适合數值運算,成為用計算機分析離散信号和系統的有力工具。對連續信号和系統,可以通過時域采樣,應用DFT進行近似譜分析。
1、用DFT對連續信号進行譜分析
1.1 實際應用DFT所面臨的問題
工程實際中,經常遇到連續信号 x a ( t ) x_a(t) xa(t),其頻譜函數 X a ( j Ω ) X_a(jΩ) Xa(jΩ)也是連續函數。為了利用DFT對 x a ( t ) x_a(t) xa(t)進行頻譜分析,先對 x a ( t ) x_a(t) xa(t)進行時域采樣,得到 x ( n ) = x a ( n T ) x(n)=x_a(nT) x(n)=xa(nT),再對 x ( n ) x(n) x(n)進行DFT,得到的 X ( k ) X(k) X(k)則是 x ( n ) x(n) x(n)的傅裡葉變換 X ( e j ω ) X(e^{j\omega}) X(ejω)在頻率區間 [ 0 , 2 π ] [0,2\pi] [0,2π]上的 N N N點等間隔采樣。這裡 x ( n ) x(n) x(n)和 X ( k ) X(k) X(k)均為有限長序列。
然而,由傅裡葉變換理論知道,若信号持續時間有限長,則其頻譜無限寬;若信号的頻譜有限寬,則其持續時間必然為無限長。是以嚴格地講,持續時間有限的帶限信号是不存在的。是以,按采樣定理采樣時,上述兩種情況下的采樣序列 x ( n ) = x a ( n T ) x(n)=x_a(nT) x(n)=xa(nT)均應為無限長,不滿足DFT的變換條件。
實際上對頻譜很寬的信号,為防止時域采樣後産生頻譜混疊失真,可用預濾波器濾除幅度較小的高頻成分,使連續信号的帶寬小于折疊頻率。對于持續時間很長的信号,采樣點數太多,以緻無法存儲和計算,隻好截取有限點進行DFT。
由上述可見,用DFT對連續信号進行頻譜分析必然是近似的,其近似程度與信号帶寬、采樣頻率和截取長度有關。實際上從工程角度看,濾除幅度很小的高頻成分和截取幅度很小的部分時間信号是允許的。是以,在下面分析中,假設 x a ( t ) x_a(t) xa(t)是經過預濾波和截取處理的有限長帶限信号。
1.2 對實際連續信号的預處理及假設
設連續信号 x a ( t ) x_a(t) xa(t)持續時間為 T p T_p Tp,最高頻率為 f c f_c fc。 x a ( t ) x_a(t) xa(t)的傅裡葉變換為 X a ( j Ω ) X_a(jΩ) Xa(jΩ),對 x a ( t ) x_a(t) xa(t)進行時域采樣得到 x ( n ) = x a ( n T ) x(n)=x_a(nT) x(n)=xa(nT), x ( n ) x(n) x(n)的傅裡葉變換為 X ( e j ω ) X(e^{j\omega}) X(ejω)。
由假設條件可知 x ( n ) x(n) x(n)的長度為
N = T p T = T p F s N=\frac{T_p}{T}=T_p F_s N=TTp=TpFs
式中, T T T為采樣間隔, F s = 1 T F_s=\frac{1}{T} Fs=T1為采樣頻率。
用 X ( k ) X(k) X(k)表示 x ( n ) x(n) x(n)的N點DFT,下面推導出 X ( k ) X(k) X(k)與 X a ( j Ω ) X_a(jΩ) Xa(jΩ)的關系,最後由此關系歸納出用 X ( k ) X(k) X(k)表示 X a ( j Ω ) X_a(jΩ) Xa(jΩ)的方法,即用DFT對連續信号進行譜分析的方法。
1.3 用DFT表示連續信号傅裡葉變換的推導過程
時域離散信号的傅裡葉變換和模拟信号傅裡葉變換之間的關系為:
X ( e j ω ) = X ^ ( j Ω ) ∣ Ω = ω T = 1 T ∑ k = − ∞ ∞ X a ( j ω − 2 π k T ) X(e^{j\omega})=\hat X(j\Omega)\Big|_{\Omega=\frac{\omega}{T}}=\frac{1}{T}\sum_{k=-∞}^{∞}X_a(j\frac{\omega-2\pi k}{T}) X(ejω)=X^(jΩ)∣∣∣Ω=Tω=T1k=−∞∑∞Xa(jTω−2πk)
由上式知道, x ( n ) x(n) x(n)的傅裡葉變換 X ( e j ω ) X(e^{j\omega}) X(ejω)與 x a ( t ) x_a(t) xa(t)的傅裡葉變換 X a ( j Ω ) X_a(jΩ) Xa(jΩ)滿足如下關系:
X ( e j ω ) = 1 T ∑ m = − ∞ ∞ X a [ j ( ω T − 2 π T m ) ] X(e^{j\omega})=\frac{1}{T}\sum_{m=-∞}^{∞}X_a\Big[j({\omega \over T}-\frac{2\pi}{T}m)\Big] X(ejω)=T1m=−∞∑∞Xa[j(Tω−T2πm)]
将 ω = Ω T \omega=\Omega T ω=ΩT代入上式,得到:
X ( e j Ω T ) = 1 T ∑ m = − ∞ ∞ X a [ j ( Ω − 2 π T m ) ] = 1 T X ~ a ( j Ω ) (3) X(e^{j\Omega T})=\frac{1}{T}\sum_{m=-∞}^{∞}X_a\Big[j(\Omega-\frac{2\pi}{T}m)\Big]=\frac{1}{T}\tilde X_a(j\Omega) \tag{3} X(ejΩT)=T1m=−∞∑∞Xa[j(Ω−T2πm)]=T1X~a(jΩ)(3)
式中
X ~ a ( j Ω ) = ∑ m = − ∞ ∞ X a [ j ( Ω − 2 π T m ) ] \tilde X_a(j\Omega)=\sum_{m=-∞}^{∞}X_a\Big[j(\Omega-\frac{2\pi}{T}m)\Big] X~a(jΩ)=m=−∞∑∞Xa[j(Ω−T2πm)]
表示模拟信号頻譜 X a ( j Ω ) X_a(j\Omega) Xa(jΩ)的周期延拓函數。
由 x ( n ) x(n) x(n)的N點DFT的定義有
X ( k ) = D F T [ x ( n ) ] N = X ( e j ω ) ∣ ω = 2 π N k k = 0 , 1 , ⋯ , N − 1 (4) X(k)=DFT[x(n)]_N=X(e^{j\omega})\Big|_{\omega=\frac{2\pi}{N}k} \qquad k=0,1,\cdots,N-1 \tag{4} X(k)=DFT[x(n)]N=X(ejω)∣∣∣ω=N2πkk=0,1,⋯,N−1(4)
将(4)式代入(3)式,得到:
X ( k ) = X ( e j 2 π N k ) = 1 T X ~ a ( j 2 π N T k ) = 1 T X ~ a ( j 2 π T p k ) k = 0 , 1 , ⋯ , N − 1 (5) X(k)=X(e^{j\frac{2\pi}{N}k})=\frac{1}{T}\tilde X_a(j\frac{2\pi}{NT}k)=\frac{1}{T}\tilde X_a(j\frac{2\pi}{T_p}k) \qquad k=0,1,\cdots,N-1 \tag{5} X(k)=X(ejN2πk)=T1X~a(jNT2πk)=T1X~a(jTp2πk)k=0,1,⋯,N−1(5)
(5)式說明了 X ( k ) X(k) X(k)與 X a ( j Ω ) X_a(jΩ) Xa(jΩ)的關系。
1.4 将以自變量ω表示的關系轉換為以頻率f為自變量
為了符合一般的頻譜描述習慣,以頻率 f f f為自變量,整理(5)式。
令
{ X a ′ ( f ) = X a ( j Ω ) ∣ Ω = 2 π f = X a ( j 2 π f ) X ~ a ′ ( f ) = X ~ a ( Ω ) ∣ Ω = 2 π f = X ~ a ( 2 π f ) \begin{cases} X'_a(f)=X_a(j \Omega)\Big|_{ \Omega=2\pi f}=X_a(j2\pi f )\\ \tilde X'_a(f)=\tilde X_a(\Omega)\Big|_{ \Omega=2\pi f}=\tilde X_a(2\pi f ) \end{cases} ⎩⎪⎨⎪⎧Xa′(f)=Xa(jΩ)∣∣∣Ω=2πf=Xa(j2πf)X~a′(f)=X~a(Ω)∣∣∣Ω=2πf=X~a(2πf)
則(5)式變為
X ( k ) = 1 T X ~ a ′ ( f ) ∣ f = k T p = 1 T X ~ a ′ ( k F ) k = 0 , 1 , ⋯ , N − 1 (5) X(k)=\frac{1}{T}\tilde X'_a(f)\Big|_{f=\frac{k}{T_p}}=\frac{1}{T}\tilde X'_a(kF) \qquad k=0,1,\cdots,N-1 \tag{5} X(k)=T1X~a′(f)∣∣∣f=Tpk=T1X~a′(kF)k=0,1,⋯,N−1(5)
由此可得
X ~ a ′ ( k F ) = T X ( k ) = T ⋅ D F T [ x ( n ) ] N k = 0 , 1 , ⋯ , N − 1 (6) \tilde X'_a(kF)=TX(k)=T\cdot DFT[x(n)]_N \qquad k=0,1,\cdots,N-1 \tag{6} X~a′(kF)=TX(k)=T⋅DFT[x(n)]Nk=0,1,⋯,N−1(6)
式中,F表示對模拟信号頻譜的采樣間隔,是以稱之為頻率分辨率, T p = N T T_p=NT Tp=NT為截斷時間長度。
F = 1 T p = 1 N T = F s N F=\frac{1}{T_p}=\frac{1}{NT}=\frac{F_s}{N} F=Tp1=NT1=NFs
1.5 對推導結果的分析
(6)式表明可以通過對連續時間采樣并進行DFT再乘以T,得到模拟信号頻譜的周期延拓函數在第一個周期 [ 0 , F s ] [0,F_s] [0,Fs]上的N點等間隔采樣 X ~ a ′ ( k F ) \tilde X'_a(kF) X~a′(kF)。
如圖所示,對滿足假設的持續時間有限的帶限信号,在滿足時域采樣定理時, X ~ a ′ ( k F ) \tilde X'_a(kF) X~a′(kF)包含了模拟信号頻譜的全部資訊( k = 0 , 1 , 2 , ⋯ , N / 2 k=0,1,2,\cdots,N/2 k=0,1,2,⋯,N/2,表示正頻率頻譜采樣; k = N / 2 + 1 , N / 2 + 2 , ⋯ , N − 1 k=N/2+1,N/2+2,\cdots,N-1 k=N/2+1,N/2+2,⋯,N−1,表示負頻率頻譜采樣)。是以上述分析方法不丢失資訊,即可由 X ( k ) X(k) X(k)恢複 X a ( j Ω ) X_a(jΩ) Xa(jΩ)或 x a ( t ) x_a(t) xa(t)。
但直接由分析結果 X ( k ) X(k) X(k)看不到 X a ( j Ω ) X_a(jΩ) Xa(jΩ)的全部頻譜特性,而隻能看到N個離散采樣點的譜線,這就是所謂的栅欄效應。對實信号,其頻譜函數具有共轭對稱性,是以分析正頻率頻譜就足夠了。不存在頻譜混疊失真時,正頻率 [ 0 , F s / 2 ] [0,F_s/2] [0,Fs/2]的頻譜采樣為
X ~ a ′ ( k F ) = T X ( k ) = T ⋅ D F T [ x ( n ) ] N k = 0 , 1 , ⋯ , N / 2 \tilde X'_a(kF)=TX(k)=T\cdot DFT[x(n)]_N \qquad k=0,1,\cdots,N/2 X~a′(kF)=TX(k)=T⋅DFT[x(n)]Nk=0,1,⋯,N/2
如果 x a ( t ) x_a(t) xa(t)持續時間無限長,上述分析中要進行截斷處理,是以會産生所謂的截斷效應,進而使譜分析産生誤差。
1.6 舉例說明截斷效應
理想低通濾波器的機關沖激響應 h a ( t ) h_a(t) ha(t)及其頻響函數 H a ( f ) H_a(f) Ha(f)如下圖:
圖中 h a ( t ) = s i n [ π ( t − α ) ] π ( t − α ) α = T p 2 h_a(t)=\frac{sin[\pi(t-\alpha)]}{\pi(t-\alpha)} \qquad \alpha=\frac{T_p}{2} ha(t)=π(t−α)sin[π(t−α)]α=2Tp
現在用DFT來分析 h a ( t ) h_a(t) ha(t)的頻率響應特性。由于 h a ( t ) h_a(t) ha(t)的持續時間為無限長,是以要截取一段 T p T_p Tp,假設 T p = 8 T_p=8 Tp=8s,采樣間隔 T = 0.25 T=0.25 T=0.25s(即采樣頻率 F s F_s Fs=4Hz),采樣點數 N = T p / T = 32 N=T_p/T=32 N=Tp/T=32,頻率采樣間隔即頻率分辨率 F = 1 / T p = 0.125 F=1/T_p=0.125 F=1/Tp=0.125Hz;由于 h a ( t ) h_a(t) ha(t)為實信号,是以僅取正頻率 [ 0 , F s / 2 ] [0,F_s/2] [0,Fs/2]頻譜采樣:
H ( k F ) = T ⋅ D F T [ h ( n ) ] 0 ≤ k ≤ 16 H(kF)=T\cdot DFT[h(n)] \qquad 0\leq k \leq 16 H(kF)=T⋅DFT[h(n)]0≤k≤16
其中, h ( n ) = h a ( n T ) R 32 ( n ) h(n)=h_a(nT)R_{32}(n) h(n)=ha(nT)R32(n)。
H ( k F ) H(kF) H(kF)如上圖(c)黑點所示。由圖可見,低頻部分近似理想低通頻響特性,而高頻誤差較大,且整個頻響都有波動。這些誤差就是由于對 h a ( t ) h_a(t) ha(t)截斷所産生的,是以通常稱之為截斷效應。為減少這種截斷誤差,可适當加長 T p T_p Tp,增加采樣點數N或用窗函數處理後再進行DFT。
1.7 譜分析範圍及頻率采樣間隔
在對連續信号進行譜分析時,主要關心兩個問題,即譜分析範圍和頻率分辨率。
譜分析範圍為 [ 0 , F s / 2 ] [0,F_s/2] [0,Fs/2],直接受采樣頻率 F s F_s Fs的限制。為了不産生頻譜混疊失真,通常要求信号的最高頻率 f c ≤ F s / 2 f_c\leq F_s/2 fc≤Fs/2。頻率分辨率用頻率采樣間隔 F F F描述, F F F表示譜分析中能夠分辨的兩個頻率分量的最小間隔。顯然, F F F越小,譜分析就越接近 X a ( j f ) X_a(jf) Xa(jf),是以 F F F越小,頻率分辨率越高。
1.8 DFT對連續信号譜分析的參數選擇原則
在已知信号的最高頻率 f c f_c fc(即譜分析範圍)時,為了避免頻率混疊現象,要求采樣速率 F s F_s Fs滿足下式:
F s > 2 f c F_s>2f_c Fs>2fc
譜分辨率 F = F s / N F=F_s/N F=Fs/N,如果保持采樣點數 N N N不變,要提高頻譜分辨率(減小 F F F)就必須降低采樣頻率,采樣頻率的降低會引起譜分析範圍變窄和頻譜混疊失真。如維持 F s F_s Fs不變,為提高頻譜分辨率可以增加采樣點數 N N N,因為 N T = T p NT=T_p NT=Tp, T = 1 / F s T=1/F_s T=1/Fs,隻有增加對信号的觀察時間 T p T_p Tp,才能增加 N N N。
T p T_p Tp和 N N N的選擇:
N > 2 f c F N>\frac{2f_c}{F} N>F2fc
T p ≥ 1 F T_p\ge \frac{1}{F} Tp≥F1
2、用DFT對離散信号(序列)進行譜分析
機關圓上的Z變換就是序列的傅裡葉變換,即
X ( e j w ) = X ( z ) ∣ z = e j w X(e^{jw})=X(z)\Big|_{z=e^{jw}} X(ejw)=X(z)∣∣∣z=ejw
X ( e j w ) X(e^{jw}) X(ejw)是 w w w的連續周期函數。如果對序列 x ( n ) x(n) x(n)進行 N N N點DFT得到 X ( k ) X(k) X(k),則 X ( k ) X(k) X(k)是在區間 [ 0 , 2 π ] [0,2\pi] [0,2π]上對 X ( e j w ) X(e^{jw}) X(ejw)的 N N N點等間隔采樣,頻譜分辨率就是采樣間隔 2 π / N 2\pi/N 2π/N。是以序列的傅裡葉變換可利用DFT(即FFT)來計算。
對周期為N的周期序列 x ~ ( n ) \tilde x(n) x~(n),根據周期序列的傅裡葉變換表達式知,其頻譜函數為
X ( e j ω ) = F T [ x ~ ( n ) ] = 2 π N ∑ k = − ∞ ∞ X ~ ( k ) δ ( ω − 2 π N k ) X(e^{j\omega})=FT[\tilde x(n)]=\frac{2\pi}{N}\sum_{k=-∞}^{∞}\tilde X(k)\delta (\omega-\frac{2\pi}{N}k) X(ejω)=FT[x~(n)]=N2πk=−∞∑∞X~(k)δ(ω−N2πk)
其中
X ~ ( k ) = D F S [ x ~ ( n ) ] = ∑ n = 0 N − 1 x ~ ( n ) e − j 2 π N k n \tilde X(k)=DFS[\tilde x(n)]=\sum_{n=0}^{N-1}\tilde x(n)e^{-j\frac{2\pi}{N}kn} X~(k)=DFS[x~(n)]=n=0∑N−1x~(n)e−jN2πkn
由于 X ~ ( k ) \tilde X(k) X~(k)以N為周期,因而 X ( e j ω ) X(e^{j\omega}) X(ejω)也是以 2 π 2\pi 2π為周期的離散譜,每個周期有N條譜線,第k條譜線位于 ω = ( 2 π / N ) k \omega=(2\pi/N)k ω=(2π/N)k處,代表 x ~ ( n ) \tilde x(n) x~(n)的k次諧波分量。而且,譜線的相對大小與 X ~ ( k ) \tilde X(k) X~(k)成正比。
由此可見,周期序列的頻譜結構可用其離散傅裡葉級數系數 X ~ ( k ) \tilde X(k) X~(k)表示。由DFT的隐含周期性知道,截取 x ~ ( n ) \tilde x(n) x~(n)的主值序列 x ( n ) = x ~ ( n ) R N ( n ) x(n)=\tilde x(n)R_N(n) x(n)=x~(n)RN(n),并進行N點DFT,得到:
X ( k ) = D F T [ x ( n ) ] N = D F T [ x ~ ( n ) R N ( n ) ] = X ~ ( k ) R N ( k ) X(k)=DFT[x(n)]_N=DFT[ \tilde x(n)R_N(n) ]=\tilde X(k)R_N(k) X(k)=DFT[x(n)]N=DFT[x~(n)RN(n)]=X~(k)RN(k)
是以可用 X ( k ) X(k) X(k)表示 x ~ ( n ) \tilde x(n) x~(n)的頻譜結構。
3、用DFT進行譜分析的誤差問題
DFT(實際中用FFT計算)可用來對連續信号和數字信号進行譜分析。在實際分析過程中,要對連續信号采樣和截斷,有些非時限資料序列也要截斷,由此可能引起分析誤差。
3.1 混疊現象
- 對連續信号進行譜分析時,首先要對其采樣,變成時域離散信号後才能用DFT(FFT)進行譜分析。
- 采樣速率 F s F_s Fs必須滿足采樣定理,否則會在 w = π w=\pi w=π(對應模拟信号 f = F s / 2 f=F_s/2 f=Fs/2)附近産生頻譜混疊現象。這時用DFT分析的結果必然在 f = F s / 2 f=F_s/2 f=Fs/2附近産生較大誤差。是以,理論上必須滿足 F s ≥ 2 f c F_s≥2f_c Fs≥2fc( f c f_c fc為連續信号的最高頻率)。
- 對 F s F_s Fs确定的情況,一般在采樣前進行預濾波,濾除高于折疊頻率 F s / 2 F_s/2 Fs/2的頻率成分,以免發生頻率混疊現象。
3.2 栅欄效應
N點DFT是在頻率區間 [ 0 , 2 π ] [0,2\pi] [0,2π]上對時域離散信号的頻譜進行N點等間隔采樣,而采樣點之間的頻譜是看不到的。這就好像從N個栅欄縫隙中觀看信号的頻譜情況,僅得到N個縫隙中看到的頻譜函數值。是以稱這種現象為栅欄效應。
由于栅欄效應,有可能漏掉(擋住)大的頻譜分量。為了把原來被“栅欄”擋住的頻譜分量檢測出來,就必須提高頻率分辨率。
(1)對有限長序列,可以在原序列尾部補零。
(2)對無限長序列,可以增大截取長度及DFT變換區間長度,進而使頻域采樣間隔變小,增加頻域采樣點數和采樣點位置,使原來漏掉的某些頻譜分量被檢測出來。
(3)對連續信号的譜分析,隻要采樣率 F s F_s Fs足夠高,且采樣點數滿足頻率分辨率要求(前述參數選擇原則),就可以認為DFT後得到的離散譜的包絡近似代表原信号的頻譜。
3.3 截斷效應
實際中遇到的序列 x ( n ) x(n) x(n)可能是無限長的,用DFT對其進行譜分析時,必須将其截短,形成有限長序列 y ( n ) = x ( n ) w ( n ) y(n)=x(n)w(n) y(n)=x(n)w(n), w ( n ) w(n) w(n)稱為窗函數,長度為N。 w ( n ) = R N ( n ) w(n)=R_N(n) w(n)=RN(n),稱為矩形窗函數。
根據傅裡葉變換的頻域卷積原理,有
Y ( e j w ) = F T [ y ( n ) ] = 1 2 π X ( e j w ) ∗ W ( e j w ) = 1 2 π ∫ − π π X ( e j θ ) W ( e j ( w − θ ) ) d θ Y(e^{jw})=FT[y(n)]=\frac{1}{2\pi}X(e^{jw})*W(e^{jw})=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(e^{j\theta})W(e^{j(w-\theta)})\,d\theta Y(ejw)=FT[y(n)]=2π1X(ejw)∗W(ejw)=2π1∫−ππX(ejθ)W(ej(w−θ))dθ
其中 X ( e j w ) = F T [ x ( n ) ] W ( e j w ) = F T [ w ( n ) ] X(e^{jw})=FT[x(n)] \qquad W(e^{jw})=FT[w(n)] X(ejw)=FT[x(n)]W(ejw)=FT[w(n)]
對矩形窗函數 w ( n ) = R N ( n ) w(n)=R_N(n) w(n)=RN(n),有
W ( e j w ) = F T [ w ( n ) ] = e − j w N − 1 2 s i n ( w N / 2 ) s i n ( w / 2 ) = W g ( w ) e j φ ( w ) W(e^{jw})=FT[w(n)]=e^{-jw\frac{N-1}{2}}\frac{sin(wN/2)}{sin(w/2)}=W_g(w)e^{j\varphi(w)} W(ejw)=FT[w(n)]=e−jw2N−1sin(w/2)sin(wN/2)=Wg(w)ejφ(w)
W g ( w ) = s i n ( w N / 2 ) s i n ( w / 2 ) W_g(w)=\frac{sin(wN/2)}{sin(w/2)} Wg(w)=sin(w/2)sin(wN/2)
例如, x ( n ) = c o s ( π 4 n ) x(n)=cos(\frac{\pi}{4}n) x(n)=cos(4πn),其頻譜為
X ( e j w ) = π ∑ l = − ∞ ∞ [ δ ( w − π 4 − 2 π l ) + δ ( w + π 4 − 2 π l ) ] X(e^{jw})=\pi \sum_{l=-∞}^{∞}\Big[\delta(w-\frac{\pi}{4}-2\pi l) +\delta(w+\frac{\pi}{4}-2\pi l) \Big] X(ejw)=πl=−∞∑∞[δ(w−4π−2πl)+δ(w+4π−2πl)]
将 x ( n ) x(n) x(n)截斷後, y ( n ) = x ( n ) R N ( n ) y(n)=x(n)R_N(n) y(n)=x(n)RN(n)的頻譜為在 δ \delta δ譜線處 W g ( w ) W_g(w) Wg(w)圖像的展寬。
由此可見,截斷後序列的頻譜 Y ( e j w ) Y(e^{jw}) Y(ejw)與原序列 X ( e j w ) X(e^{jw}) X(ejw)必然有差别,這種差别對譜分析的影響主要表現在如下兩個方面:
3.3.1 洩露
原來序列 x ( n ) x(n) x(n)的頻譜是離散譜線,經截斷後,使原來的離散譜線向附近展寬,通常稱這種展寬為洩露。顯然,洩露使頻譜變模糊,使譜分辨率降低。
從上圖中可以看出,頻譜洩露程度與窗函數幅度譜的主瓣寬度直接相關,而在所有的窗函數中,矩形窗的主瓣是最窄的,但其旁瓣的幅度也最大。是以,在窗函數長度N相同時,用矩形窗截取,産生的洩露最小。
3.3.2 譜間幹擾
在主譜線兩邊形成很多旁瓣,引起不同頻率分量間的幹擾(簡稱譜間幹擾),特别是強信号譜的旁瓣可能湮沒弱信号的主譜線,或者把強信号譜的旁瓣誤以為是另一頻率的信号的譜線,進而造成假信号,這樣就會使譜分析産生較大偏差。由于矩形窗的旁瓣最大,是以,用矩形窗截取時,産生的譜間幹擾最大。
洩露與譜間幹擾
由于洩露和譜間幹擾是由信号截斷引起的,是以稱之為截斷效應。矩形窗 ∣ ω ∣ < 2 π N |\omega|<\frac{2\pi}{N} ∣ω∣<N2π的部分為主瓣,增加N可使 W g ( ω ) W_g(\omega) Wg(ω)的主瓣變窄,減小洩露,提高頻率分辨率,但旁瓣的相對幅度并不減小。為了減小譜間幹擾,應用其它形狀的窗函數 w ( n ) w(n) w(n)代替矩形窗。但在N一定時,旁瓣幅度越小的窗函數,其主瓣就越寬。是以,在DFT變換區間(即截取長度)N一定時,隻能以降低譜分析分辨率為代價,換取譜間幹擾的減小。
栅欄效應與頻率分辨率
栅欄效應與頻率分辨率是不同的兩個概念。如果截取長度為N的一段資料序列,則可以在其後面補N個零,再進行2N點DFT,使栅欄寬度減半,進而減輕了栅欄效應。
但是,這種截短後補零的方法不能提高頻率分辨率。因為截短已經使頻譜變模糊,補零後僅使采樣間隔變小,但得到的頻譜采樣仍是已經變模糊的頻譜,是以頻率分辨率沒有提高。是以,要提高頻率分辨率,就必須對原始信号截取的長度加長(對模拟信号,就是增加采樣時間 T p T_p Tp的長度)。