維格納分布 (又名韋格納分佈 ,英文: Wigner Distribution Function ,縮寫為WDF ) 是由1963年的諾貝爾物理學獎 得主尤金·维格纳 ,于1932年首次引用的一個新的方程式 。
眾所皆知,傅立葉變換 對於研究穩態(時間獨立)的訊號(波形)是一項非常有用的工具,然而,訊號(波形)一般來說在時間上並非是獨立的,這樣的訊號或是波形傅立葉變換並無法有效地完全分析其特性,因此對於一個非穩態的訊號完全分析需要測量出時間以及頻率上的表現。本頁面介紹的數學函數是時頻分析 中的基礎方法,在1980年,Claasen,Mecklenbrauker對WDF做了更進一步的研究。除此之外,線性時頻分析中,STFT、Gabor transform和WDF扮演了相當重要的角色,其中WDF對於分析很多非穩態的隨機訊號都有很好的表現,例如:量子力學 、光學 、聲學 、通訊 、生物工程 、訊號處理 和影像處理 。有時也被用在分析地震的資料,以及處理聲音的相位失真。
定義 維格納分布有許多不同的定義,而此處的定義是特別針對時頻分析而定的。若給定一時間序列 x [ t ] {\displaystyle x[t]} ,它的非平穩自相關函數 如下公式所列
C x ( t 1 , t 2 ) = ⟨ ( x [ t 1 ] − μ [ t 1 ] ) ( x [ t 2 ] − μ [ t 2 ] ) ∗ ⟩ , {\displaystyle C_{x}(t_{1},t_{2})=\left\langle \left(x[t_{1}]-\mu [t_{1}]\right)\left(x[t_{2}]-\mu [t_{2}]\right)^{*}\right\rangle ,}
其中 ⟨ ⋯ ⟩ {\displaystyle \langle \cdots \rangle } 代表所有可能實驗的程序的平均, μ [ t ] {\displaystyle \mu [t]} 代表平均,其可能是時間的函數也有可能不是。維格納函數 W x ( t , f ) {\displaystyle W_{x}(t,f)} 起初是以包含時間平均 t = ( t 1 + t 2 ) / 2 {\displaystyle t=(t_{1}+t_{2})/2} 與時間差 τ = t 1 − t 2 {\displaystyle \tau =t_{1}-t_{2}} 的自相關函數和時間差進行傅立葉轉換來表示,如下:
W x ( t , f ) = ∫ − ∞ ∞ C x ( t + τ 2 , t − τ 2 ) e − 2 π i τ f d τ . {\displaystyle W_{x}(t,f)=\int _{-\infty }^{\infty }C_{x}\left(t+{\frac {\tau }{2}},t-{\frac {\tau }{2}}\right)\,e^{-2\pi i\tau f}\,d\tau .}
對於單一零平均的時間序列,維格納函數可以簡化如下:
定義一 W x ( t , f ) = ∫ − ∞ ∞ x ( t + τ 2 ) x ∗ ( t − τ 2 ) e − j 2 π τ f d τ {\displaystyle W_{x}(t,f)=\int _{-\infty }^{\infty }x(t+{\frac {\tau }{2}})x^{*}(t-{\frac {\tau }{2}})e^{-j2\pi \tau f}\,d\tau } .....(1)
定義二 W x ( t , ω ) = ∫ − ∞ ∞ x ( t + τ 2 ) x ∗ ( t − τ 2 ) e − j ω τ d τ {\displaystyle W_{x}(t,\omega )=\int _{-\infty }^{\infty }x(t+{\frac {\tau }{2}})x^{*}(t-{\frac {\tau }{2}})e^{-j\omega \tau }\,d\tau } .....(2)定義二與定義一之間的關係 : ω = 2 π f {\displaystyle \omega =2\pi f}
其他定義 1 2 π W x ( t , ω ) = ∫ − ∞ ∞ x ( t + τ 2 ) x ∗ ( t − τ 2 ) e − j ω τ d τ {\displaystyle {\sqrt {\frac {1}{2\pi }}}W_{x}(t,\omega )=\int _{-\infty }^{\infty }x(t+{\frac {\tau }{2}})x^{*}(t-{\frac {\tau }{2}})e^{-j\omega \tau }\,d\tau } .....(3)在聲納 和雷達 系統中,傳送出去的聲波的反射波可以用來偵測目標物的位置跟速度,在很多情形下,收到的訊號因為都普勒位移,所以跟原本的訊號並不一樣。Woodward(1953) 改寫了原本的公式
A x ( t , ω ) = ∫ − ∞ ∞ x ( τ + t 2 ) x ∗ ( τ − t 2 ) e − j ω τ d τ {\displaystyle A_{x}(t,\omega )=\int _{-\infty }^{\infty }x(\tau +{\frac {t}{2}})x^{*}(\tau -{\frac {t}{2}})e^{-j\omega \tau }\,d\tau } 這個公式被稱為Woodward ambiguity function,這個式子在雷達系統的訊號處理和設計上扮演重要的角色。
而維格納分布亦為科恩系列分布 的其中一種特例,當科恩系列分布中的 Φ ( η , τ ) = 1 {\displaystyle \Phi (\eta ,\tau )=1} 時,科恩系列分布會是維格納分布。
WDF和STFT的比較 WDF、STFT和Gabor transform 都佔了時頻分析中非常重要的地位,在這邊比較一下它們之間的差別。
WDF STFT 清晰度 較好 較差 相交項的問題 嚴重 無 複雜度 高 低 處理隨機程序 可 不可
相交項其實就是處理的過程中產生的額外訊號,是不想要的,WDF的清晰度和複雜度是彼此做取捨的,可以依不同的情況或是不同的方法來決定是否要使用WDF或是另外兩種。
WDF的優缺點 在這裡列出WDF主要的優缺點
優點 :
1.有良好的解析度,尤其是對單一成分,且瞬時頻率變化不為2次式以上。
2.有好的數學運算性質(見WDF的數學性質)。
3.可用於分析隨機程序(見WDF與隨機程序的關係)。
缺點 :
1.有相交項(cross term)的問題,改進方法請見 改進型韋格納分佈 。
2.需要更多的時間去計算,若訊號時間越長,則需要更久的時間。
3.不是一對一函數,無法辨別相位部分,例如: W D F [ x ( t ) ] = W D F [ x ( t ) e j ϕ ] {\displaystyle WDF[x(t)]=WDF[x(t){e^{j\phi }}]}
4.不適合分析瞬時頻率變化為2次式以上的型態,即 e j t n , n ≠ 0 , 1 , 2 {\displaystyle {e^{j{t^{n}}}},n\neq 0,1,2} 。
相交項特性
WDF與隨機程序的關係 對於一個隨機程序x(t),我們無法得知其確切的值,因此會將其值表示為一個機率函數,通常E[x(t)] = 0 for any t
將x(t)的維格納分布取期望值後可得其譜密度(Power spectral density,PSD),如下公式所列
E [ W x ( t , f ) ] = ∫ − ∞ ∞ E [ x ( t + τ / 2 ) x ∗ ( t − τ / 2 ) ] ⋅ e − j 2 π f τ ⋅ d τ {\displaystyle E[W_{x}(t,f)]=\textstyle \int \limits _{-\infty }^{\infty }\displaystyle E[x(t+\tau /2)x^{*}(t-\tau /2)]\cdot e^{-j2\pi f\tau }\cdot d\tau }
= ∫ − ∞ ∞ R x ( t , τ ) ⋅ e − j 2 π f τ ⋅ d τ = S x ( t , f ) {\displaystyle =\textstyle \int \limits _{-\infty }^{\infty }\displaystyle R_{x}(t,\tau )\cdot e^{-j2\pi f\tau }\cdot d\tau =S_{x}(t,f)}
當x(t)的統計特性不隨時間變化時,可稱x(t)為平穩的隨機程序,其譜密度也可簡化為 E [ W x ( t , f ) ] = S x ( f ) {\displaystyle E[W_{x}(t,f)]=S_{x}(f)} ,也就是說維格納函數能初略的告訴我們譜密度如何隨時間進行變化。維格納函數能在平穩程序對所有時間t都簡化成譜密度,然而也等同於非平穩的自相關函數,這也是維格納分布的動機。
下圖為一個平穩的隨機程序進行維格納分析後的例子,可明顯看出此信號不隨時間變化,也就是時頻分析結果為水平線。反之,亦可利用時頻分析結果是否為水平線判斷該訊號是否為一平穩的訊號。
而在訊號處理中常見的白雜訊 ,其譜密度 S x ( f ) = σ {\displaystyle S_{x}(f)=\sigma } ,其中 σ {\displaystyle \sigma } 為一個常數。白雜訊的維格納分布如下圖,可看出此雜訊在所有時間及頻率都存在著。
維格納分布的相交項在處理隨機程序時派上用場,相對的,沒有相交項的短時距傅立葉轉換,則無法用於隨機程序,如下公式所示,只有在零平均隨機程序時, E [ X ( t , f ) ] = 0 {\displaystyle E[X(t,f)]=0}
E [ X ( t , f ) ] = E [ ∫ t − B t + B x ( τ ) w ( t − τ ) e − j 2 π f τ d τ ] = ∫ t − B t + B E [ x ( τ ) ] w ( t − τ ) e − j 2 π f τ d τ {\displaystyle E[X(t,f)]=E[\textstyle \int \limits _{t-B}^{t+B}\displaystyle x(\tau )w(t-\tau )e^{-j2\pi f\tau }d\tau ]=\textstyle \int \limits _{t-B}^{t+B}\displaystyle E[x(\tau )]w(t-\tau )e^{-j2\pi f\tau }d\tau }
常見的時頻分析例子 以下的例子說明如何用WDF來做時頻分析
常數訊號 輸入訊號為常數,則時頻分佈為一條線重合於時軸,如果'x( t) = 1,則:
W x ( t , f ) = ∫ − ∞ ∞ e − i 2 π τ f d τ = δ ( f ) . {\displaystyle W_{x}(t,f)=\int _{-\infty }^{\infty }e^{-i2\pi \tau \,f}\,d\tau =\delta (f).}
弦波訊號 輸入訊號為弦波,則時頻分佈為一條線平行於時軸,如果x (t ) = e i2πkt ,則:
W x ( t , f ) = ∫ − ∞ ∞ e i 2 π k ( t + τ 2 ) e − i 2 π k ( t − τ 2 ) e − i 2 π τ f d τ = ∫ − ∞ ∞ e − i 2 π τ ( f − k ) d τ = δ ( f − k ) . {\displaystyle {\begin{aligned}W_{x}(t,f)&=\int _{-\infty }^{\infty }e^{i2\pi k\left(t+{\frac {\tau }{2}}\right)}e^{-i2\pi k\left(t-{\frac {\tau }{2}}\right)}e^{-i2\pi \tau \,f}\,d\tau \\&=\int _{-\infty }^{\infty }e^{-i2\pi \tau \left(f-k\right)}\,d\tau \\&=\delta (f-k).\end{aligned}}}
啁啾聲信號 啁啾聲訊號的瞬時頻率隨時間線性,表示時頻分佈為一條斜值線,例如
x ( t ) = e i 2 π k t 2 {\displaystyle x(t)=e^{i2\pi kt^{2}}} ,則瞬時頻率為:
1 2 π d ( 2 π k t 2 ) d t = 2 k t , {\displaystyle {\frac {1}{2\pi }}{\frac {d(2\pi kt^{2})}{dt}}=2kt~,} 故WDF為:
W x ( t , f ) = ∫ − ∞ ∞ e i 2 π k ( t + τ 2 ) 2 e − i 2 π k ( t − τ 2 ) 2 e − i 2 π τ f d τ = ∫ − ∞ ∞ e i 4 π k t τ e − i 2 π τ f d τ = ∫ − ∞ ∞ e − i 2 π τ ( f − 2 k t ) d τ = δ ( f − 2 k t ) . {\displaystyle {\begin{aligned}W_{x}(t,f)&=\int _{-\infty }^{\infty }e^{i2\pi k\left(t+{\frac {\tau }{2}}\right)^{2}}e^{-i2\pi k\left(t-{\frac {\tau }{2}}\right)^{2}}e^{-i2\pi \tau \,f}\,d\tau \\&=\int _{-\infty }^{\infty }e^{i4\pi kt\tau }e^{-i2\pi \tau f}\,d\tau \\&=\int _{-\infty }^{\infty }e^{-i2\pi \tau (f-2kt)}\,d\tau \\&=\delta (f-2kt)~.\end{aligned}}}
余弦訊號 x(t) = cos(440 π {\displaystyle \pi } t), 当 t 小於 0.5, 頻率 f = 220Hz x(t) = cos(660 π {\displaystyle \pi } t), 当 0.5 小於等於 t 小於 1, 頻率 f = 330Hz x(t) = cos(524 π {\displaystyle \pi } t), 当 t 大於等於 1, 頻率 f = 262Hz
單位脈衝訊號 因為單位脈衝包含所有的頻率分佈,且在時間不等於零時沒值,故WDF為通過原點的且與時軸垂直的線
W x ( t , f ) = ∫ − ∞ ∞ δ ( t + τ 2 ) δ ( t − τ 2 ) e − i 2 π τ f d τ = 4 ∫ − ∞ ∞ δ ( 2 t + τ ) δ ( 2 t − τ ) e − i 2 π τ f d τ = 4 δ ( 4 t ) e i 4 π t f = δ ( t ) e i 4 π t f = δ ( t ) . {\displaystyle {\begin{aligned}W_{x}(t,f)&=\int _{-\infty }^{\infty }\delta \left(t+{\frac {\tau }{2}}\right)\delta \left(t-{\frac {\tau }{2}}\right)e^{-i2\pi \tau \,f}\,d\tau \\&=4\int _{-\infty }^{\infty }\delta (2t+\tau )\delta (2t-\tau )e^{-i2\pi \tau f}\,d\tau \\&=4\delta (4t)e^{i4\pi tf}\\&=\delta (t)e^{i4\pi tf}\\&=\delta (t).\end{aligned}}}
方波訊號 x ( t ) = { 1 | t | < 1 / 2 0 otherwise {\displaystyle x(t)={\begin{cases}1&|t|<1/2\\0&{\text{otherwise}}\end{cases}}\qquad } , W x ( t , f ) = 1 π f sin ( f [ 1 − 2 | t | ] ) {\displaystyle W_{x}(t,f)={\frac {1}{\pi f}}\sin(f[1-2|t|])} .
WDF的數學性質 (1)投射特性 | x ( t ) | 2 = ∫ − ∞ ∞ W x ( t , f ) d f {\displaystyle |x(t)|^{2}=\int _{-\infty }^{\infty }W_{x}(t,f)\,df} , | X ( f ) | 2 = ∫ − ∞ ∞ W x ( t , f ) d t {\displaystyle |X(f)|^{2}=\int _{-\infty }^{\infty }W_{x}(t,f)\,dt} (2)能量特性 | x ( t ) | 2 = ∬ − ∞ ∞ W x ( t , f ) d t d f = ∫ − ∞ ∞ | x ( t ) | 2 d t = ∫ − ∞ ∞ | X ( f ) | 2 d f {\displaystyle |x(t)|^{2}=\iint _{-\infty }^{\infty }W_{x}(t,f)\,dt\,df=\int _{-\infty }^{\infty }|x(t)|^{2}\,dt=\int _{-\infty }^{\infty }|X(f)|^{2}\,df} (3)回覆特性 ∫ − ∞ ∞ W x ( t 2 , f ) e j 2 π f t d f = x ( t ) ∙ x ∗ ( 0 ) {\displaystyle \int _{-\infty }^{\infty }W_{x}({\frac {t}{2}},f)e^{j2\pi ft}\,df=x(t)\bullet x^{*}(0)} , ∫ − ∞ ∞ W x ( t , f 2 ) e − j 2 π f t d t = X ( f ) ∙ X ∗ ( 0 ) {\displaystyle \int _{-\infty }^{\infty }W_{x}(t,{\frac {f}{2}})e^{-j2\pi ft}\,dt=X(f)\bullet X^{*}(0)} (4)Mean 條件 If x ( t ) = | x ( t ) | e j 2 π ϕ ( t ) , X ( f ) = | X ( f ) | e j 2 π Ψ ( f ) {\displaystyle x(t)=\left\vert x(t)\right\vert e^{j2\pi \phi (t)},\ X(f)=\left\vert X(f)\right\vert e^{j2\pi \Psi (f)}} then ϕ ′ ( t ) = | x ( t ) | − 2 ∫ − ∞ ∞ f × W x ( t , f ) d f {\displaystyle \phi '(t)={\left\vert x(t)\right\vert }^{-2}\textstyle \int \limits _{-\infty }^{\infty }\displaystyle f\times W_{x}(t,f)\ df} , − Ψ ′ ( f ) = | X ( f ) | − 2 ∫ − ∞ ∞ t × W x ( t , f ) d t {\displaystyle -\Psi '(f)={\left\vert X(f)\right\vert }^{-2}\textstyle \int \limits _{-\infty }^{\infty }\displaystyle t\times W_{x}(t,f)\ dt}
(5)Moment特性 ∬ − ∞ ∞ t n W x ( t , f ) d t d f = ∫ − ∞ ∞ t n | x ( t ) | 2 d t {\displaystyle \iint _{-\infty }^{\infty }t^{n}W_{x}(t,f)\,dt\,df=\int _{-\infty }^{\infty }t^{n}|x(t)|^{2}\,dt} , ∬ − ∞ ∞ f n W x ( t , f ) d t d f = ∫ − ∞ ∞ f n | X ( f ) | 2 d t {\displaystyle \iint _{-\infty }^{\infty }f^{n}W_{x}(t,f)\,dt\,df=\int _{-\infty }^{\infty }f^{n}|X(f)|^{2}\,dt} (6) W x ( t , f ) {\displaystyle W_{x}(t,f)} 是實數 W x ( t , f ) ¯ = W x ( t , f ) {\displaystyle {\overline {W_{x}(t,f)}}=W_{x}(t,f)} (7)區域特性 If x ( t ) = 0 {\displaystyle x(t)=0} for t > t 0 {\displaystyle t>t_{0}} then W x ( t , f ) = 0 {\displaystyle W_{x}(t,f)=0} for t > t 0 {\displaystyle t>t_{0}} , If x ( t ) = 0 {\displaystyle x(t)=0} for t < t 0 {\displaystyle t<t_{0}} then W x ( t , f ) = 0 {\displaystyle W_{x}(t,f)=0} for t < t 0 {\displaystyle t<t_{0}} (8)乘法特性 If y ( t ) = x ( t ) h ( t ) {\displaystyle y(t)=x(t)h(t)} ,then W y ( t , f ) = ∫ − ∞ ∞ W x ( t , ρ ) W h ( t − ρ , f ) d ρ {\displaystyle W_{y}(t,f)=\int _{-\infty }^{\infty }W_{x}(t,\rho )W_{h}(t-\rho ,f)\,d\rho } (9)摺積特性 If y ( t ) = ∫ − ∞ ∞ x ( t − τ ) h ( τ ) d τ {\displaystyle y(t)=\int _{-\infty }^{\infty }x(t-\tau )h(\tau )\,d\tau } ,then W y ( t , f ) = ∫ − ∞ ∞ W x ( ρ , f ) W h ( t − ρ , f ) d ρ {\displaystyle W_{y}(t,f)=\int _{-\infty }^{\infty }W_{x}(\rho ,f)W_{h}(t-\rho ,f)\,d\rho } (10)相關特性 If y ( t ) = ∫ − ∞ ∞ x ( t + τ ) h ∗ ( τ ) d τ {\displaystyle y(t)=\int _{-\infty }^{\infty }x(t+\tau )h^{*}(\tau )\,d\tau } ,then W y ( t , ω ) = ∫ − ∞ ∞ W x ( ρ , ω ) W h ( − t + ρ , ω ) d ρ {\displaystyle W_{y}(t,\omega )=\int _{-\infty }^{\infty }W_{x}(\rho ,\omega )W_{h}(-t+\rho ,\omega )\,d\rho } (11)時間平移特性 If y ( t ) = x ( t − t 0 ) {\displaystyle y(t)=x(t-t_{0})} , then W y ( t , f ) = W x ( t − t 0 , f ) {\displaystyle W_{y}(t,f)=W_{x}(t-t_{0},f)} (12)調變特性 If y ( t ) = e j 2 π f 0 t x ( t ) {\displaystyle y(t)=e^{j2\pi f_{0}t}x(t)} , then W y ( t , f ) = W x ( t , f − f 0 ) {\displaystyle W_{y}(t,f)=W_{x}(t,f-f_{0})}
WDF的數學性質證明
WDF實現方法 以下為電腦計算WDF的實現方式
直接運算(暴力法) 複雜度: T F ( 2 Q + 1 ) {\displaystyle TF(2Q+1)} 使用離散傅立葉變換 複雜度: T N log 2 N {\displaystyle TN{\log _{2}}N} 使用Chirp-Z 轉換 複雜度 : T N log 2 N {\displaystyle TN{\log _{2}}N} ,通常為使用離散傅立葉變換 的2~3倍,但限制比使用離散傅立葉變換 少 在使用這三個方法前,先來做個前提討論
從定義一出發
W x ( t , f ) = ∫ − ∞ ∞ x ( t + τ / 2 ) ⋅ x ∗ ( t − τ / 2 ) e − j 2 π τ f ⋅ d τ {\displaystyle {W_{x}}\left({t,f}\right)=\int _{-\infty }^{\infty }{x\left({t+\tau /2}\right)\cdot }{x^{*}}\left({t-\tau /2}\right)\,{e^{-j2\pi \,\tau \,f}}\cdot d\tau }
令 τ ′ = τ / 2 {\displaystyle \tau '=\tau /2}
W x ( t , f ) = 2 ∫ − ∞ ∞ x ( t + τ ′ ) ⋅ x ∗ ( t − τ ′ ) e − j 4 π τ ′ f ⋅ d τ ′ {\displaystyle {W_{x}}\left({t,f}\right)=2\int _{-\infty }^{\infty }{x\left({t+\tau '}\right)\cdot }{x^{*}}\left({t-\tau '}\right)\,{e^{-j4\pi \,\tau '\,f}}\cdot d\tau '}
再令 t = n Δ t , f = m Δ f , τ ′ = p Δ t {\displaystyle t=n\Delta _{t},f=m\Delta _{f},\tau '=p\Delta _{t}} ,則上述式子則為
W x ( n Δ t , m Δ f ) = 2 ∑ p = − ∞ ∞ x ( ( n + p ) Δ t ) x ∗ ( ( n − p ) Δ t ) exp ( − j 4 π m p Δ t Δ f ) Δ t {\displaystyle {W_{x}}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=2\sum \limits _{p=-\infty }^{\infty }{x\left({(n+p){\Delta _{t}}}\right){x^{*}}\left({(n-p){\Delta _{t}}}\right)\exp \left({-j4\pi \,mp{\Delta _{t}}{\Delta _{f}}}\right){\Delta _{t}}}}
下面介紹的三種方法都是從這條式子開始推導
注意事項 :
若x(t)是無限長的訊號,則p要從負無限加到正無限,這點不易實現。
若x(t)為有限長的訊號,則p範圍可以縮小,就可能實現。
故下面三種方法都是在第2種情況下討論,即x(t)為有限長訊號,p範圍可以縮小
我們假設 x ( t ) = 0 f o r t < n 1 Δ t a n d t > n 2 Δ t {\displaystyle \ x(t)=0\ \ for\ \ t<n_{1}\Delta _{t}\ \ and\ \ t>n_{2}\Delta _{t}}
直接運算(暴力法) 限制條件 :
只有一個 : 要滿足Nyquist criterion
Δ t < 1 2 B {\displaystyle {\Delta _{t}}<{\frac {1}{2B}}} ,其中B是 x ( t + τ ) x ∗ ( t − τ ) {\displaystyle x\left({t+\tau }\right){x^{*}}\left({t-\tau }\right)\,} 的頻寬,大約是x(t)的兩倍。
推導 :
x ( t ) = 0 f o r t < n 1 Δ t a n d t > n 2 Δ t {\displaystyle \ x(t)=0\ \ for\ \ t<n_{1}\Delta _{t}\ \ and\ \ t>n_{2}\Delta _{t}}
所以當 n + p ∉ [ n 1 , n 2 ] o r n − p ∉ [ n 1 , n 2 ] {\displaystyle n+p\notin [{n_{1}},{n_{2}}]{\rm {orn-p}}\notin {\rm {[}}{n_{1}},{n_{2}}{\rm {]}}} 時,
x ( ( n + p ) Δ t ) x ∗ ( ( n − p ) Δ t ) = 0 {\displaystyle x\left({(n+p){\Delta _{t}}}\right){x^{*}}\left({(n-p){\Delta _{t}}}\right)=0}
固定中間的n值 ( n Δ t {\displaystyle n\Delta _{t}} ) 來探討p的範圍
n 1 ≤ n + p ≤ n 2 → n 1 − n ≤ p ≤ n 2 − n {\displaystyle {n_{1}}\leq n+p\leq {n_{2}}\to {n_{1}}-n\leq p\leq {n_{2}}-n}
即 max ( n 1 − n , n − n 2 ) ≤ p ≤ min ( n 2 − n , n − n 1 ) {\displaystyle \max({n_{1}}-n,n-{n_{2}})\leq p\leq \min({n_{2}}-n,n-{n_{1}})} -– (1)
n 1 ≤ n − p ≤ n 2 → n 1 − n ≤ − p ≤ n 2 − n → n − n 2 ≤ p ≤ n − n 1 {\displaystyle {n_{1}}\leq n-p\leq {n_{2}}\to {\rm {}}{{\rm {n}}_{1}}-n\leq -p\leq {n_{2}}-n\to n-{n_{2}}\leq p\leq n-{n_{1}}}
即 − min ( n 2 − n , n − n 1 ) ≤ p ≤ min ( n 2 − n , n − n 1 ) {\displaystyle -\min({n_{2}}-n,n-{n_{1}})\leq p\leq \min({n_{2}}-n,n-{n_{1}})} -- (2)
其中 (1) & (2) 的下限是同義的
故(1) & (2)皆可改寫為
− min ( n 2 − n , n − n 1 ) ≤ p ≤ min ( n 2 − n , n − n 1 ) {\displaystyle -\min({n_{2}}-n,n-{n_{1}})\leq p\leq \min({n_{2}}-n,n-{n_{1}})}
且可以發現 ( n 2 − n ) Δ t , ( n − n 1 ) Δ t {\displaystyle ({n_{2}}-n){\Delta _{t}},(n-{n_{1}}){\Delta _{t}}} 代表 n Δ t {\displaystyle n\Delta _{t}} 離兩個邊界的距離
注意事項: 當 n > n2 或 n < n1 時,將沒有 p 能滿足上面的不等式
最後推導出的式子如下
W x ( n Δ t , m Δ f ) = 2 ∑ p = − Q Q x ( ( n + p ) Δ t ) x ∗ ( ( n − p ) Δ t ) exp ( − j 4 π m p Δ t Δ f ) Δ t {\displaystyle {W_{x}}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=2\sum \limits _{p=-Q}^{Q}{x\left({(n+p){\Delta _{t}}}\right){x^{*}}\left({(n-p){\Delta _{t}}}\right)\exp \left({-j4\pi \,mp{\Delta _{t}}{\Delta _{f}}}\right){\Delta _{t}}}}
其中 Q = min ( n 2 − n , n − n 1 ) , p ∈ [ − Q , Q ] , n ∈ [ n 1 , n 2 ] {\displaystyle Q=\min({n_{2}}-n,n-{n_{1}})\,,p\in [-Q,Q],n\in [{n_{1}},{n_{2}}]}
使用離散傅立葉變換 限制條件 :
(1)要滿足Nyquist criterion
Δ t < 1 2 B {\displaystyle {\Delta _{t}}<{\frac {1}{2B}}} ,其中B是 x ( t + τ ) x ∗ ( t − τ ) {\displaystyle x\left({t+\tau }\right){x^{*}}\left({t-\tau }\right)\,} 的頻寬,大約是x(t)的兩倍。
(2) Δ t Δ f = 1 2 N {\displaystyle {\Delta _{t}}{\Delta _{f}}={\textstyle {1 \over {2N}}}}
(3) N ≥ 2 Q + 1 {\displaystyle N\geq 2Q+1}
推導 :
前提討論的式子可以改寫為
W x ( n Δ t , m Δ f ) = 2 Δ t ∑ p = − Q Q x ( ( n + p ) Δ t ) x ∗ ( ( n − p ) Δ t ) e − j 2 π m p N {\displaystyle {W_{x}}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=2{\Delta _{t}}\sum \limits _{p=-Q}^{Q}{x\left({(n+p){\Delta _{t}}}\right){x^{*}}\left({(n-p){\Delta _{t}}}\right){e^{-j{\textstyle {{2\pi \,mp} \over N}}}}}}
令 q = p + Q → p = q − Q {\displaystyle q=p+Q\to p=q-Q}
W x ( n Δ t , m Δ f ) = 2 Δ t e j 2 π m Q N ∑ q = 0 2 Q x ( ( n + q − Q ) Δ t ) x ∗ ( ( n − q + Q ) Δ t ) e − j 2 π m q N {\displaystyle {W_{x}}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=2{\Delta _{t}}{e^{j{\textstyle {{2\pi \,mQ} \over N}}}}\sum \limits _{q=0}^{2Q}{x\left({(n+q-Q){\Delta _{t}}}\right){x^{*}}\left({(n-q+Q){\Delta _{t}}}\right){e^{-j{\textstyle {{2\pi \,mq} \over N}}}}}}
針對中間 x ( ( n + q − Q ) Δ t ) x ∗ ( ( n − q + Q ) Δ t ) {\displaystyle x\left({(n+q-Q){\Delta _{t}}}\right){x^{*}}\left({(n-q+Q){\Delta _{t}}}\right)} 項
令
c 1 ( q ) = x ( ( n + q − Q ) Δ t ) x ∗ ( ( n − q + Q ) Δ t ) , f o r 0 ≤ q ≤ 2 Q {\displaystyle {c_{1}}\left(q\right)=x\left({(n+q-Q){\Delta _{t}}}\right){x^{*}}\left({(n-q+Q){\Delta _{t}}}\right){\rm {}},for{\rm {0}}\leq q\leq 2{\rm {Q}}}
c 1 ( q ) = 0 , f o r 2 Q + 1 ≤ q ≤ N − 1 {\displaystyle {c_{1}}\left(q\right)=0{\rm {}}\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad ,for{\rm {2}}Q+1\leq q\leq N-1}
最後得出的式子如下
W x ( n Δ t , m Δ f ) = 2 Δ t e j 2 π m Q N ∑ q = 0 N − 1 c 1 ( q ) e − j 2 π m q N {\displaystyle {W_{x}}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=2{\Delta _{t}}{e^{j{\textstyle {{2\pi \,mQ} \over N}}}}\sum \limits _{q=0}^{N-1}{{c_{1}}\left(q\right){e^{-j{\textstyle {{2\pi \,mq} \over N}}}}}}
其中
Q = min ( n 2 − n , n − n 1 ) , n ∈ [ n 1 , n 2 ] {\displaystyle Q=\min({n_{2}}-n,n-{n_{1}})\,,n\in [{n_{1}},{n_{2}}]}
c 1 ( q ) = x ( ( n + q − Q ) Δ t ) x ∗ ( ( n − q + Q ) Δ t ) , f o r 0 ≤ q ≤ 2 Q {\displaystyle {c_{1}}\left(q\right)=x\left({(n+q-Q){\Delta _{t}}}\right){x^{*}}\left({(n-q+Q){\Delta _{t}}}\right){\rm {}},for{\rm {0}}\leq q\leq 2{\rm {Q}}}
c 1 ( q ) = 0 , f o r 2 Q + 1 ≤ q ≤ N − 1 {\displaystyle {c_{1}}\left(q\right)=0{\rm {}}\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad ,for{\rm {2}}Q+1\leq q\leq N-1}
使用Chirp-Z 轉換 限制條件 :
只有一個 : 要滿足Nyquist criterion
Δ t < 1 2 B {\displaystyle {\Delta _{t}}<{\frac {1}{2B}}} ,其中B是 x ( t + τ ) x ∗ ( t − τ ) {\displaystyle x\left({t+\tau }\right){x^{*}}\left({t-\tau }\right)\,} 的頻寬,大約是x(t)的兩倍。
推導 :
前提討論的式子可改寫為
W x ( n Δ t , m Δ f ) = 2 Δ t e − j 2 π m 2 Δ t Δ f ∑ p = − Q Q x ( ( n + p ) Δ t ) x ∗ ( ( n − p ) Δ t ) e − j 2 π p 2 Δ t Δ f e j 2 π ( p − m ) 2 Δ t Δ f {\displaystyle {W_{x}}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=2{\Delta _{t}}\,{e^{-j2\pi \,{m^{2}}{\Delta _{t}}{\Delta _{f}}}}\sum \limits _{p=-Q}^{Q}{x\left({(n+p){\Delta _{t}}}\right){x^{*}}\left({(n-p){\Delta _{t}}}\right){e^{-j2\pi \,{p^{2}}{\Delta _{t}}{\Delta _{f}}}}{e^{j2\pi \,{{(p-m)}^{2}}{\Delta _{t}}{\Delta _{f}}}}}}
計算分成3步驟
STEP 1 : x 1 ( n , p ) = x ( ( n + p ) Δ t ) x ∗ ( ( n − p ) Δ t ) e − j 2 π p 2 Δ t Δ f {\displaystyle {x_{1}}\left({n,p}\right)=x\left({(n+p){\Delta _{t}}}\right){x^{*}}\left({(n-p){\Delta _{t}}}\right){e^{-j2\pi \,{p^{2}}{\Delta _{t}}{\Delta _{f}}}}}
STEP 2 : X 2 [ n , m ] = ∑ p = n − Q n + Q x 1 [ p ] c [ m − p ] {\displaystyle {X_{2}}\left[{n,m}\right]=\sum \limits _{p=n-Q}^{n+Q}{{x_{1}}\left[p\right]\,c\left[{m-p}\right]}} , 其中 c [ m ] = e j 2 π m 2 Δ t Δ f {\displaystyle c\left[m\right]={e^{j2\pi \,{m^{2}}{\Delta _{t}}{\Delta _{f}}}}}
STEP 3 : X ( n Δ t , m Δ f ) = 2 Δ t e − j 2 π m 2 Δ t Δ f X 2 [ n , m ] {\displaystyle X\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=2{\Delta _{t}}\,{e^{-j2\pi \;{m^{2}}{\Delta _{t}}{\Delta _{f}}}}{X_{2}}\left[{n,m}\right]}
延伸變化
視窗型韋格納分佈 視窗型韋格納分佈 (Windowed Wigner Distribution Function ),在韋格納分佈 中,當x(t)為無限長訊號時,WDF很難去實現它。所以在積分中加入一個新的函數 ,目的是擷取x(t)中的片段來計算,不需從負無限積分到正無限。
定義 W x ( t , f ) = ∫ − ∞ ∞ w ( τ ) x ( t + τ / 2 ) ⋅ x ∗ ( t − τ / 2 ) e − j 2 π τ f ⋅ d τ {\displaystyle {W_{x}}\left({t,f}\right)=\int _{-\infty }^{\infty }{w\left(\tau \right)x\left({t+\tau /2}\right)\cdot }{x^{*}}\left({t-\tau /2}\right)\,{e^{-j2\pi \,\tau \,f}}\cdot d\tau } , 其中 w ( τ ) {\displaystyle w(\tau )} 為實數且為有限長訊號
原始韋格納分佈 定義 W x ( t , f ) = ∫ − ∞ ∞ x ( t + τ / 2 ) ⋅ x ∗ ( t − τ / 2 ) e − j 2 π τ f ⋅ d τ {\displaystyle {W_{x}}\left({t,f}\right)=\int _{-\infty }^{\infty }{x\left({t+\tau /2}\right)\cdot }{x^{*}}\left({t-\tau /2}\right)\,{e^{-j2\pi \,\tau \,f}}\cdot d\tau }
優缺點 降低運算時間,因為 w ( τ ) {\displaystyle w(\tau )} 為有限長函數。 可以有效降低相交項(cross term)問題,但不能完全消除(詳見下方說明)。 一些相交項(cross term)問題仍被保留。 可能不符合譜密度 (Power spectral density)的定義。 一些好用的數學運算性質會消失。
實現方法 從定義出發 W x ( t , f ) = ∫ − ∞ ∞ w ( τ ) x ( t + τ / 2 ) ⋅ x ∗ ( t − τ / 2 ) e − j 2 π τ f ⋅ d τ {\displaystyle {W_{x}}\left({t,f}\right)=\int _{-\infty }^{\infty }{w\left(\tau \right)x\left({t+\tau /2}\right)\cdot }{x^{*}}\left({t-\tau /2}\right)\,{e^{-j2\pi \,\tau \,f}}\cdot d\tau } 令 τ = τ ′ / 2 {\displaystyle \tau =\tau '/2} W x ( t , f ) = 2 ∫ − ∞ ∞ w ( 2 τ ′ ) x ( t + τ ′ ) ⋅ x ∗ ( t − τ ′ ) e − j 4 π τ ′ f ⋅ d τ ′ {\displaystyle {W_{x}}\left({t,f}\right)=2\int _{-\infty }^{\infty }{w\left({2\tau '}\right)x\left({t+\tau '}\right)\cdot }{x^{*}}\left({t-\tau '}\right)\,{e^{-j4\pi \,\tau '\,f}}\cdot d\tau '} 再令 t = n Δ t , f = m Δ f , τ ′ = p Δ t {\displaystyle t=n\Delta _{t},f=m\Delta _{f},\tau '=p\Delta _{t}} W x ( n Δ t , m Δ f ) = 2 ∑ p = − ∞ ∞ w ( 2 p Δ t ) x ( ( n + p ) Δ t ) x ∗ ( ( n − p ) Δ t ) e − j 4 π m p Δ t Δ f Δ t {\displaystyle {W_{x}}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=2\sum \limits _{p=-\infty }^{\infty }{w\left({2p{\Delta _{t}}}\right)x\left({(n+p){\Delta _{t}}}\right){x^{*}}\left({(n-p){\Delta _{t}}}\right){e^{-j4\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}}} 假設w(t) = 0 for |t| > B 即 w ( 2 p Δ t ) = 0 f o r p < − Q ∧ p > Q {\displaystyle w\left({2p{\Delta _{t}}}\right)=0{\rm {}}\ for{\rm {}}\ p<-Q{\rm {}}\ \land \ p>Q} 其中 Q = B 2 Δ t {\displaystyle Q={\frac {B}{2{\Delta _{t}}}}} 如此一來,p範圍便可縮小。 W x ( n Δ t , m Δ f ) = 2 ∑ p = − Q Q w ( 2 p ) x ( ( n + p ) Δ t ) x ∗ ( ( n − p ) Δ t ) e − j 4 π m p Δ t Δ f Δ t {\displaystyle {W_{x}}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=2\sum \limits _{p=-Q}^{Q}{w\left({2p}\right)x\left({(n+p){\Delta _{t}}}\right){x^{*}}\left({(n-p){\Delta _{t}}}\right){e^{-j4\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}}}
避免相交項的原因 從定義出發 W x ( t , f ) = ∫ − ∞ ∞ w ( τ ) x ( t + τ / 2 ) ⋅ x ∗ ( t − τ / 2 ) e − j 2 π τ f ⋅ d τ {\displaystyle {W_{x}}\left({t,f}\right)=\int _{-\infty }^{\infty }{w\left(\tau \right)x\left({t+\tau /2}\right)\cdot }{x^{*}}\left({t-\tau /2}\right)\,{e^{-j2\pi \,\tau \,f}}\cdot d\tau } ,其中 w ( τ ) {\displaystyle w(\tau )} 為實數且為有限長訊號
假設 x ( t ) = δ ( t − t 1 ) + δ ( t − t 2 ) {\displaystyle x(t){\rm {=}}\delta (t-{t_{1}})+\delta (t-{t_{2}})} 的情況下,比較有無mask function所產生的不同結果
x(t)示意圖 理想情形 : W x ( t , f ) = 0 f o r t ≠ t 1 , t 2 {\displaystyle {W_{x}}(t,f){\rm {=0}}\ for{\rm {}}\ t\neq {t_{1}},{t_{2}}}
沒有使用mask function 即mask function w ( τ ) = 1 {\displaystyle w(\tau )=1}
W x ( t , f ) = ∫ − ∞ ∞ x ( t + τ / 2 ) ⋅ x ∗ ( t − τ / 2 ) e − j 2 π τ f ⋅ d τ = ∫ − ∞ ∞ [ δ ( t + τ 2 − t 1 ) + δ ( t + τ 2 − t 2 ) ] ⋅ [ δ ( t − τ 2 − t 1 ) + δ ( t − τ 2 − t 2 ) ] e − j 2 π τ f ⋅ d τ = 4 ∫ − ∞ ∞ [ δ ( τ + 2 t − 2 t 1 ) + δ ( τ + 2 t − 2 t 2 ) ] ⏟ f i r s t t e r m ⋅ [ δ ( τ − 2 t + 2 t 1 ) + δ ( τ − 2 t + 2 t 2 ) ] ⏟ s e c o n d t e r m e − j 2 π τ f ⋅ d τ {\displaystyle {\begin{array}{l}{W_{x}}\left({t,f}\right)=\int _{-\infty }^{\infty }{x\left({t+\tau /2}\right)\cdot }{x^{*}}\left({t-\tau /2}\right)\,{e^{-j2\pi \,\tau \,f}}\cdot d\tau \\=\int _{-\infty }^{\infty }{\left[{\delta \left({t+{\frac {\tau }{2}}-{t_{1}}}\right)+\delta \left({t+{\frac {\tau }{2}}-{t_{2}}}\right)}\right]\cdot }\left[{\delta \left({t-{\frac {\tau }{2}}-{t_{1}}}\right)+\delta \left({t-{\frac {\tau }{2}}-{t_{2}}}\right)}\right]{e^{-j2\pi \,\tau \,f}}\cdot d\tau \\=4\int _{-\infty }^{\infty }\underbrace {\left[{\delta \left({\tau +2t-2{t_{1}}}\right)+\delta \left({\tau +2t-2{t_{2}}}\right)}\right]} _{first\ term}\cdot \underbrace {\left[{\delta \left({\tau -2t+2{t_{1}}}\right)+\delta \left({\tau -2t+2{t_{2}}}\right)}\right]} _{second\ term}{e^{-j2\pi \,\tau \,f}}\cdot d\tau \end{array}}}
Wx(t,f)示意圖 總共有3種情況要討論,如下圖,可見cross term在沒有使用mask function時,無法被消除
ideal x(t)。Auto term 為自相關項。 Cross term為相交項
使用mask function W x ( t , f ) = ∫ − ∞ ∞ w ( τ ) x ( t + τ / 2 ) ⋅ x ∗ ( t − τ / 2 ) e − j 2 π τ f ⋅ d τ = 4 ∫ − ∞ ∞ w ( τ ) [ δ ( τ + 2 t − 2 t 1 ) + δ ( τ + 2 t − 2 t 2 ) ] ⏟ f i r s t t e r m ⋅ [ δ ( τ − 2 t + 2 t 1 ) + δ ( τ − 2 t + 2 t 2 ) ] ⏟ s e c o n d t e r m e − j 2 π τ f ⋅ d τ {\displaystyle {\begin{array}{l}{W_{x}}\left({t,f}\right)=\int _{-\infty }^{\infty }w(\tau ){x\left({t+\tau /2}\right)\cdot }{x^{*}}\left({t-\tau /2}\right)\,{e^{-j2\pi \,\tau \,f}}\cdot d\tau \\=4\int _{-\infty }^{\infty }w(\tau )\underbrace {\left[{\delta \left({\tau +2t-2{t_{1}}}\right)+\delta \left({\tau +2t-2{t_{2}}}\right)}\right]} _{first\ term}\cdot \underbrace {\left[{\delta \left({\tau -2t+2{t_{1}}}\right)+\delta \left({\tau -2t+2{t_{2}}}\right)}\right]} _{second\ term}{e^{-j2\pi \,\tau \,f}}\cdot d\tau \end{array}}}
假設 w ( τ ) = 0 f o r | τ | > B > 0 {\displaystyle w(\tau )=0{\rm {}}\ for|\tau |>B>0} ,且 B < t 2 − t 1 {\displaystyle B{\rm {<}}{{\rm {t}}_{2}}{\rm {-}}{{\rm {t}}_{1}}}
由於 w ( τ ) {\displaystyle w(\tau )} 只在-B到B有值,故乘上 w ( τ ) {\displaystyle w(\tau )} 就能去除相交項(Cross term),只保留下圖中兩條紅線中間的區域,也就是Auto terms。
ideal x(t)。 但上述其實是理想的情況,x(t)為窄頻信號Delta function
如果X(τ)寬度太寬或是有ripple的話,Cross term仍會有殘留,示意圖如下
non ideal x(t)。 藍色線為X(τ)的訊號,若X(τ)的寬度太寬或是有ripple產生,就可能會跑進 w ( τ ) {\displaystyle w(\tau )} 的範圍裡面,進而導致無法完全濾除Cross term。
總結 cross term 只有在訊號每個成分的寬度都小於2B,且時間差 t 2 − t 1 {\displaystyle t_{2}-t_{1}} 都大於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.