天天看點

信号處理--離散傅裡葉變換(DFT)

一、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−1​x(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−1​x(n)WNkn​k=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)]=N1​k=0∑N−1​X(k)WN−kn​n=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−1​x(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−1​x(n)WNkn​k=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π​k​k=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π​k​k=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)
%%%%%%%%%% 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是一種時域和頻域均離散化的變換,适合數值運算,成為用計算機分析離散信号和系統的有力工具。對連續信号和系統,可以通過時域采樣,應用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​​=Tp​Fs​

式中, 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ω​​=T1​k=−∞∑∞​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ω)=T1​m=−∞∑∞​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)=T1​m=−∞∑∞​Xa​[j(Ω−T2π​m)]=T1​X~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π​k​k=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)=T1​X~a​(jNT2π​k)=T1​X~a​(jTp​2π​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)=T1​X~a′​(f)∣∣∣​f=Tp​k​​=T1​X~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)]N​k=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=Tp​1​=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)。

信号處理--離散傅裡葉變換(DFT)

如圖所示,對滿足假設的持續時間有限的帶限信号,在滿足時域采樣定理時, 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)]N​k=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)

現在用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−1​x~(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π1​X(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−1​sin(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)​

信号處理--離散傅裡葉變換(DFT)

例如, 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)圖像的展寬。

信号處理--離散傅裡葉變換(DFT)

由此可見,截斷後序列的頻譜 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​的長度)。

繼續閱讀