維格納分佈

来自testwiki
跳转到导航 跳转到搜索

維格納分布(又名韋格納分佈,英文: Wigner Distribution Function,縮寫為WDF) 是由1963年的諾貝爾物理學獎得主尤金·维格纳,于1932年首次引用的一個新的方程式

眾所皆知,傅立葉變換對於研究穩態(時間獨立)的訊號(波形)是一項非常有用的工具,然而,訊號(波形)一般來說在時間上並非是獨立的,這樣的訊號或是波形傅立葉變換並無法有效地完全分析其特性,因此對於一個非穩態的訊號完全分析需要測量出時間以及頻率上的表現。本頁面介紹的數學函數是時頻分析中的基礎方法,在1980年,Claasen,Mecklenbrauker對WDF做了更進一步的研究。除此之外,線性時頻分析中,STFT、Gabor transform和WDF扮演了相當重要的角色,其中WDF對於分析很多非穩態的隨機訊號都有很好的表現,例如:量子力學光學聲學通訊生物工程訊號處理影像處理。有時也被用在分析地震的資料,以及處理聲音的相位失真。

定義

維格納分布有許多不同的定義,而此處的定義是特別針對時頻分析而定的。若給定一時間序列x[t],它的非平穩自相關函數如下公式所列

   Cx(t1,t2)=(x[t1]μ[t1])(x[t2]μ[t2])*,

其中代表所有可能實驗的程序的平均,μ[t]代表平均,其可能是時間的函數也有可能不是。維格納函數Wx(t,f)起初是以包含時間平均t=(t1+t2)/2與時間差τ=t1t2的自相關函數和時間差進行傅立葉轉換來表示,如下:

  Wx(t,f)=Cx(t+τ2,tτ2)e2πiτfdτ.

對於單一零平均的時間序列,維格納函數可以簡化如下:

定義一

Wx(t,f)=x(t+τ2)x*(tτ2)ej2πτfdτ.....(1)

定義二

Wx(t,ω)=x(t+τ2)x*(tτ2)ejωτdτ.....(2)

定義二與定義一之間的關係 : ω=2πf

其他定義

12πWx(t,ω)=x(t+τ2)x*(tτ2)ejωτdτ.....(3)

聲納雷達系統中,傳送出去的聲波的反射波可以用來偵測目標物的位置跟速度,在很多情形下,收到的訊號因為都普勒位移,所以跟原本的訊號並不一樣。Woodward(1953) 改寫了原本的公式

Ax(t,ω)=x(τ+t2)x*(τt2)ejωτdτ

這個公式被稱為Woodward ambiguity function,這個式子在雷達系統的訊號處理和設計上扮演重要的角色。

而維格納分布亦為科恩系列分布的其中一種特例,當科恩系列分布中的Φ(η,τ)=1時,科恩系列分布會是維格納分布。

WDF、STFT和Gabor的比較

WDF、STFT和Gabor都佔了時頻分析中非常重要的地位,在這邊比較一下WDF、STFT、Gabor和WDF的變體Cohen's Class Distribution之間的差別。

WDF STFT Gabor Cohen
清晰度 較好 較差 適中 偏差
相交項的問題 嚴重 輕微 嚴重
複雜度
處理隨機程序 不可
頻率分辨率
時間分辨率
能量集中性 較好 較差 適中 適中
數值穩定性 較差 較好 較好 較差

清晰度:描述時頻圖在時間和頻率域上解析信號細節的能力。

相交項的問題:多分量信號時可能出現不想要的相交項(cross term),導致時頻圖難以解讀。

複雜度:指運算複雜度,較高的複雜度意味著需要更多的計算資源。

處理隨機程序:對具有隨機性質的信號進行有效分析和表徵的能力。

頻率分辨率:頻域解析度的能力,越高越能分辨接近的頻率分量。

時間分辨率:時間域解析度的能力,越高越能定位信號的時間變化。

能量集中性:時頻圖中能量分布是否聚集在主要分量,分散則會影響分析效果。

數值穩定性:算法在不同條件下(如隨機噪聲或非穩態信號)保持穩定表現的能力。

WDF的優缺點

在這裡列出WDF主要的優缺點

優點 :

1.有良好的解析度,尤其是對單一成分,且瞬時頻率變化不為2次式以上,適合分析瞬態與多變信號。

2.有良好的數學運算性質(見WDF的數學性質)。

3.對隨機過程和非穩態信號表現良好,可用於分析隨機程序(見WDF與隨機程序的關係)。

4.信號能量分佈高度集中於真實分量,便於解析信號結構。

5.在時間和頻率域上的邊際分佈與原始信號一致, 滿足投影性質。

6.提供準確的能量分佈估計,便於信號重建。

7.能準確表示信號瞬時頻率,有助於動態信號分析,有限制(缺點4)。

8.無窗函數限制,能解析接近頻率分量。

9.對非線性信號具有優勢,便於非線性行為檢測。

10.對含有噪聲或隨機成分的信號也能表現出色。

缺點 :

1.有相交項(cross term)的問題,降低時頻圖可解讀性,改進方法請見 改進型韋格納分佈

2.需執行二維積分或卷積,對計算資源要求高。若訊號時間越長,則需要更久的時間。

3.不是一對一函數,無法辨別相位部分,例如: WDF[x(t)]=WDF[x(t)ejϕ]

4.不適合分析瞬時頻率變化為2次式以上的型態,即ejtn,n0,1,2

5.低信噪比環境下,噪聲可能被放大,影響準確性。

6.不滿足線性運算,影響某些應用的便利性。

7.結果分佈可能對輸入條件敏感,缺乏唯一性。

相交項特性

WDF不是一個線性的轉換,由於x(t)的signal auto-correlation function x(t+τ2)x*(tτ2),如果有兩種以上不同性質的訊號疊加,會產生相交項。然而相交項卻有重要且有用的物理意義,像是可以用來分析期望值,相反的,短時距傅立葉轉換就沒有此特性,詳見維格納分佈與隨機程序的關係。以下數學方程式對於WDF後會產生相交項。

  • x(t)={cos(2πt)t2cos(4πt)2<t2cos(3πt)t>2


  • x(t)=eit3

WDF與隨機程序的關係

對於一個隨機程序x(t),我們無法得知其確切的值,因此會將其值表示為一個機率函數,通常E[x(t)] = 0 for any t

將x(t)的維格納分布取期望值後可得其譜密度(Power spectral density,PSD),如下公式所列

E[Wx(t,f)]=E[x(t+τ/2)x*(tτ/2)]ej2πfτdτ

=Rx(t,τ)ej2πfτdτ=Sx(t,f)

當x(t)的統計特性不隨時間變化時,可稱x(t)為平穩的隨機程序,其譜密度也可簡化為E[Wx(t,f)]=Sx(f),也就是說維格納函數能初略的告訴我們譜密度如何隨時間進行變化。維格納函數能在平穩程序對所有時間t都簡化成譜密度,然而也等同於非平穩的自相關函數,這也是維格納分布的動機。

下圖為一個平穩的隨機程序進行維格納分析後的例子,可明顯看出此信號不隨時間變化,也就是時頻分析結果為水平線。反之,亦可利用時頻分析結果是否為水平線判斷該訊號是否為一平穩的訊號。

而在訊號處理中常見的白雜訊,其譜密度Sx(f)=σ,其中σ為一個常數。白雜訊的維格納分布如下圖,可看出此雜訊在所有時間及頻率都存在著。

維格納分布的相交項在處理隨機程序時派上用場,相對的,沒有相交項的短時距傅立葉轉換,則無法用於隨機程序,如下公式所示,只有在零平均隨機程序時,E[X(t,f)]=0

E[X(t,f)]=E[tBt+Bx(τ)w(tτ)ej2πfτdτ]=tBt+BE[x(τ)]w(tτ)ej2πfτdτ

常見的時頻分析例子

以下的例子說明如何用WDF來做時頻分析

常數訊號

輸入訊號為常數,則時頻分佈為一條線重合於時軸,如果'x(t) = 1,則:

Wx(t,f)=ei2πτfdτ=δ(f).

弦波訊號

輸入訊號為弦波,則時頻分佈為一條線平行於時軸,如果Template:Math,則:

Wx(t,f)=ei2πk(t+τ2)ei2πk(tτ2)ei2πτfdτ=ei2πτ(fk)dτ=δ(fk).

啁啾聲信號

啁啾聲訊號的瞬時頻率隨時間線性,表示時頻分佈為一條斜值線,例如

x(t)=ei2πkt2 ,

則瞬時頻率為:

12πd(2πkt2)dt=2kt,

故WDF為:

Wx(t,f)=ei2πk(t+τ2)2ei2πk(tτ2)2ei2πτfdτ=ei4πktτei2πτfdτ=ei2πτ(f2kt)dτ=δ(f2kt).

余弦訊號

x(t) = cos(440πt), 当 t 小於 0.5, 頻率 f = 220Hz
x(t) = cos(660πt), 当 0.5 小於等於 t 小於 1, 頻率 f = 330Hz
x(t) = cos(524πt), 当 t 大於等於 1, 頻率 f = 262Hz

單位脈衝訊號

因為單位脈衝包含所有的頻率分佈,且在時間不等於零時沒值,故WDF為通過原點的且與時軸垂直的線

Wx(t,f)=δ(t+τ2)δ(tτ2)ei2πτfdτ=4δ(2t+τ)δ(2tτ)ei2πτfdτ=4δ(4t)ei4πtf=δ(t)ei4πtf=δ(t).


方波訊號

x(t)={1|t|<1/20otherwise ,
Wx(t,f)=1πfsin(f[12|t|]) .

WDF的數學性質

(1)投射特性 |x(t)|2=Wx(t,f)df, |X(f)|2=Wx(t,f)dt
(2)能量特性 |x(t)|2=Wx(t,f)dtdf=|x(t)|2dt=|X(f)|2df
(3)回覆特性 Wx(t2,f)ej2πftdf=x(t)x*(0), Wx(t,f2)ej2πftdt=X(f)X*(0)
(4)Mean 條件 If x(t)=|x(t)|ej2πϕ(t), X(f)=|X(f)|ej2πΨ(f)

then ϕ(t)=|x(t)|2f×Wx(t,f) df, Ψ(f)=|X(f)|2t×Wx(t,f) dt

(5)Moment特性 tnWx(t,f)dtdf=tn|x(t)|2dt, fnWx(t,f)dtdf=fn|X(f)|2dt
(6)Wx(t,f)是實數 Wx(t,f)=Wx(t,f)
(7)區域特性 If x(t)=0 for t>t0 then Wx(t,f)=0 for t>t0, If x(t)=0 for t<t0 then Wx(t,f)=0 for t<t0
(8)乘法特性 If y(t)=x(t)h(t),then Wy(t,f)=Wx(t,ρ)Wh(tρ,f)dρ
(9)摺積特性 If y(t)=x(tτ)h(τ)dτ,then Wy(t,f)=Wx(ρ,f)Wh(tρ,f)dρ
(10)相關特性 If y(t)=x(t+τ)h*(τ)dτ,then Wy(t,ω)=Wx(ρ,ω)Wh(t+ρ,ω)dρ
(11)時間平移特性 If y(t)=x(tt0), then Wy(t,f)=Wx(tt0,f)
(12)調變特性 If y(t)=ej2πf0tx(t), then Wy(t,f)=Wx(t,ff0)

WDF的數學性質證明

WDF滿足永遠是實數的性質,以下是證明:

Wx(t,f)=x(t+τ2)x*(tτ2)ej2πfτdτ
Wx(t,f)=x*(t+τ2)x(tτ2)ej2πfτdτ
τ=τ代入,變數乘上負號,因此積分範圍會變成
Wx(t,f)=x*(tτ2)x(t+τ2)ej2πfτ(dτ)=x*(tτ2)x(t+τ2)ej2πfτdτ=x(t+τ2)x*(tτ2)ej2πfτ=Wx(t,f)
故WDF永遠是實數

WDF實現方法

以下為電腦計算WDF的實現方式

  1. 直接運算(暴力法) 複雜度:TF(2Q+1)
  2. 使用離散傅立葉變換 複雜度:TNlog2N
  3. 使用Chirp-Z 轉換 複雜度 : TNlog2N,通常為使用離散傅立葉變換的2~3倍,但限制比使用離散傅立葉變換

在使用這三個方法前,先來做個前提討論

從定義一出發

Wx(t,f)=x(t+τ/2)x*(tτ/2)ej2πτfdτ

τ=τ/2

Wx(t,f)=2x(t+τ)x*(tτ)ej4πτfdτ

再令 t=nΔt,f=mΔf,τ=pΔt ,則上述式子則為

Wx(nΔt,mΔf)=2p=x((n+p)Δt)x*((np)Δt)exp(j4πmpΔtΔf)Δt

下面介紹的三種方法都是從這條式子開始推導

注意事項 :

若x(t)是無限長的訊號,則p要從負無限加到正無限,這點不易實現。

若x(t)為有限長的訊號,則p範圍可以縮小,就可能實現。


故下面三種方法都是在第2種情況下討論,即x(t)為有限長訊號,p範圍可以縮小

我們假設 x(t)=0  for  t<n1Δt  and  t>n2Δt


直接運算(暴力法)


限制條件 :

只有一個 : 要滿足Nyquist criterion

Δt<12B,其中B是x(t+τ)x*(tτ)的頻寬,大約是x(t)的兩倍。

推導 :

 x(t)=0  for  t<n1Δt  and  t>n2Δt

所以當n+p[n1,n2]ornp[n1,n2] 時,

x((n+p)Δt)x*((np)Δt)=0

固定中間的n值 (nΔt) 來探討p的範圍

n1n+pn2n1npn2n

max(n1n,nn2)pmin(n2n,nn1)-– (1)

n1npn2n1npn2nnn2pnn1

min(n2n,nn1)pmin(n2n,nn1)-- (2)

其中 (1) & (2) 的下限是同義的

故(1) & (2)皆可改寫為

min(n2n,nn1)pmin(n2n,nn1)

且可以發現 (n2n)Δt,(nn1)Δt 代表 nΔt 離兩個邊界的距離

注意事項: 當 n > n2 或 n < n1 時,將沒有 p 能滿足上面的不等式

最後推導出的式子如下

Wx(nΔt,mΔf)=2p=QQx((n+p)Δt)x*((np)Δt)exp(j4πmpΔtΔf)Δt

其中 Q=min(n2n,nn1),p[Q,Q],n[n1,n2]

使用離散傅立葉變換


限制條件 :

(1)要滿足Nyquist criterion

Δt<12B,其中B是x(t+τ)x*(tτ)的頻寬,大約是x(t)的兩倍。

(2)ΔtΔf=12N

(3) N2Q+1

推導 :

前提討論的式子可以改寫為

Wx(nΔt,mΔf)=2Δtp=QQx((n+p)Δt)x*((np)Δt)ej2πmpN

q=p+Qp=qQ

Wx(nΔt,mΔf)=2Δtej2πmQNq=02Qx((n+qQ)Δt)x*((nq+Q)Δt)ej2πmqN

針對中間x((n+qQ)Δt)x*((nq+Q)Δt)

c1(q)=x((n+qQ)Δt)x*((nq+Q)Δt),for0q2Q

c1(q)=0,for2Q+1qN1


最後得出的式子如下

Wx(nΔt,mΔf)=2Δtej2πmQNq=0N1c1(q)ej2πmqN

其中

Q=min(n2n,nn1),n[n1,n2]

c1(q)=x((n+qQ)Δt)x*((nq+Q)Δt),for0q2Q

c1(q)=0,for2Q+1qN1

使用Chirp-Z 轉換


限制條件 :

只有一個 : 要滿足Nyquist criterion

Δt<12B,其中B是x(t+τ)x*(tτ)的頻寬,大約是x(t)的兩倍。

推導 :

前提討論的式子可改寫為

Wx(nΔt,mΔf)=2Δtej2πm2ΔtΔfp=QQx((n+p)Δt)x*((np)Δt)ej2πp2ΔtΔfej2π(pm)2ΔtΔf

計算分成3步驟

STEP 1 : x1(n,p)=x((n+p)Δt)x*((np)Δt)ej2πp2ΔtΔf

STEP 2 : X2[n,m]=p=nQn+Qx1[p]c[mp] , 其中c[m]=ej2πm2ΔtΔf

STEP 3 : X(nΔt,mΔf)=2Δtej2πm2ΔtΔfX2[n,m]

延伸變化

視窗型韋格納分佈

視窗型韋格納分佈Template:Lang),在韋格納分佈中,當x(t)為無限長訊號時,WDF很難去實現它。所以在積分中加入一個新的函數 ,目的是擷取x(t)中的片段來計算,不需從負無限積分到正無限。

定義

Wx(t,f)=w(τ)x(t+τ/2)x*(tτ/2)ej2πτfdτ , 其中w(τ)為實數且為有限長訊號

原始韋格納分佈定義 Wx(t,f)=x(t+τ/2)x*(tτ/2)ej2πτfdτ

優缺點

  • 優點 :
  1. 降低運算時間,因為w(τ)為有限長函數。
  2. 可以有效降低相交項(cross term)問題,但不能完全消除(詳見下方說明)。
  • 缺點 :
  1. 一些相交項(cross term)問題仍被保留。
  2. 可能不符合譜密度(Power spectral density)的定義。
  3. 一些好用的數學運算性質會消失。

實現方法

從定義出發 Wx(t,f)=w(τ)x(t+τ/2)x*(tτ/2)ej2πτfdττ=τ/2 Wx(t,f)=2w(2τ)x(t+τ)x*(tτ)ej4πτfdτ 再令t=nΔt,f=mΔf,τ=pΔt Wx(nΔt,mΔf)=2p=w(2pΔt)x((n+p)Δt)x*((np)Δt)ej4πmpΔtΔfΔt 假設w(t) = 0 for |t| > B 即 w(2pΔt)=0 for p<Q  p>Q 其中Q=B2Δt 如此一來,p範圍便可縮小。 Wx(nΔt,mΔf)=2p=QQw(2p)x((n+p)Δt)x*((np)Δt)ej4πmpΔtΔfΔt

避免相交項的原因

從定義出發 Wx(t,f)=w(τ)x(t+τ/2)x*(tτ/2)ej2πτfdτ,其中w(τ)為實數且為有限長訊號

假設x(t)=δ(tt1)+δ(tt2)的情況下,比較有無mask function所產生的不同結果

x(t)示意圖

理想情形 : Wx(t,f)=0 for tt1,t2


沒有使用mask function

即mask function w(τ)=1

Wx(t,f)=x(t+τ/2)x*(tτ/2)ej2πτfdτ=[δ(t+τ2t1)+δ(t+τ2t2)][δ(tτ2t1)+δ(tτ2t2)]ej2πτfdτ=4[δ(τ+2t2t1)+δ(τ+2t2t2)]first term[δ(τ2t+2t1)+δ(τ2t+2t2)]second termej2πτfdτ

Wx(t,f)示意圖

總共有3種情況要討論,如下圖,可見cross term在沒有使用mask function時,無法被消除

ideal x(t)。Auto term 為自相關項。 Cross term為相交項


使用mask function

Wx(t,f)=w(τ)x(t+τ/2)x*(tτ/2)ej2πτfdτ=4w(τ)[δ(τ+2t2t1)+δ(τ+2t2t2)]first term[δ(τ2t+2t1)+δ(τ2t+2t2)]second termej2πτfdτ

假設w(τ)=0 for|τ|>B>0,且B<t2t1

由於w(τ)只在-B到B有值,故乘上w(τ)就能去除相交項(Cross term),只保留下圖中兩條紅線中間的區域,也就是Auto terms。

ideal x(t)。

但上述其實是理想的情況,x(t)為窄頻信號Delta function

如果X(τ)寬度太寬或是有ripple的話,Cross term仍會有殘留,示意圖如下

non ideal x(t)。

藍色線為X(τ)的訊號,若X(τ)的寬度太寬或是有ripple產生,就可能會跑進w(τ)的範圍裡面,進而導致無法完全濾除Cross term。

總結

cross term 只有在訊號每個成分的寬度都小於2B,且時間差t2t1都大於B時,才能被消除

此方法可以消除相交項(cross term)。

消除相交項(cross term)問題,在某些情況下比加伯轉換擁有更好的清晰度。

参见

參考書目、資料來源

  • Jian-Jiun Ding, Time frequency analysis and wavelet transform class note, the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2018.
  • Jian-Jiun Ding, Time-frequency analysis and wavelet transform class note, the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2021.
  • Jian-Jiun Ding, Time-frequency analysis and wavelet transform class note, the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2023.
  • Jian-Jiun Ding, Time-frequency analysis and wavelet transform class note, the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2024.