韋格納分佈的時頻分析

韋格納分布Wigner Distribution Function(縮寫為WDF)，是時頻分析中的一種分析方式，以下為其方程式的轉換：

$W_x(t,f)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}x(t+\tau /2)\cdot x^{*}(t-\tau /2)e^{-j2\pi f}\cdot d\tau$

經整理後得到 $W_x(t,f)=2\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}x(t+{\tau }^{\prime })\cdot x^{*}(t-{\tau }^{\prime })e^{-j4\pi {\tau }^{\prime }f}\cdot d{\tau }^{\prime },using{\tau }^{\prime }=\tau /2$

連續訊號要得到其數位實現，必須經過採樣(sampling), 使得 $t=n{\mathrm{\Delta }}_t,f=m{\mathrm{\Delta }}_f,{\tau }^{\prime }=p{\mathrm{\Delta }}_t$

重新回顧上述式子可以得到，$W_x(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=2\sum _{p=-\mathrm{\infty }}^{\mathrm{\infty }}x((n+p){\mathrm{\Delta }}_t)\cdot x^{*}((n-p){\mathrm{\Delta }}_t)e^{-j4\pi mp{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{\mathrm{\Delta }}_t$

當 $x(t)$ 不是 time-limited signal，會很難實現。

於是，通常我們會假設 $x(t)=0fort<n_1{\mathrm{\Delta }}_{t}andt>n_2{\mathrm{\Delta }}_t$，也就是 $x(t)$ 為一個有限長度的訊號：

如圖所示：

此時 $x((n+p){\mathrm{\Delta }}_t)x^{*}((n-p){\mathrm{\Delta }}_t)=0,if(n+p)\in ̸[n_1,n_2]or(n-p)\in ̸[n1,n2]$

繼續討論 $p$ 的範圍 ( 於$n$為一固定值時 )：

1. 對於$(n+p)$：$n_1\le (n+p)\le n_2⟶(n_1-n)\le p\le (n_2-n)$
2. 對於$(n-p)$：$n_1\le (n-p)\le n_2⟶(n_1-n)\le -p\le (n_2-n)$，也就是 $(n-n_2)\le p\le (n-n_1)$

於是便得知 $p$ 的範圍為：$\max (n_1-n,n-n_2)\le p\le \min (n_2-n,n-n_1)⟶-\min (n_2-n,n-n_1)\le p\le \min (n_2-n,n-n_1)⟶-Q\le p\le Q,Q=min(n_2-n,n-n_1)$

用圖片來理解：

注意：當 $n>n_2$ 或 $n<n_1$時，找不到適當的 $p$ 來滿足不等式。

整理前述說明後，WDF的數位實現的數學式可整理為：

$ifx(t)=0fort<n_1{\mathrm{\Delta }}_{t}andt>n_2{\mathrm{\Delta }}_t$，

$W_x(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=2\sum _{p=-Q}^{Q}x((n+p){\mathrm{\Delta }}_t)x^{*}((n-p){\mathrm{\Delta }}_t)e^{-j4\pi mp{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{\mathrm{\Delta }}_t$

(

$Q=\min (n_2-n,n-n_1),(varieswithn)$ $and$ $p\in [-Q,Q],n\in [n_1,n_2]$

)

以下提供3種實現方式：

1. 暴力法 (Direct Implementation)
2. 離散傅立葉轉換(Using Discrete time Fourier Transform)
3. Chirp-Z轉換

暴力法(Direct Implementation)

根據 $W_x(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=2\sum _{p=-Q}^{Q}x((n+p){\mathrm{\Delta }}_t)x^{*}((n-p){\mathrm{\Delta }}_t)e^{-j4\pi mp{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{\mathrm{\Delta }}_t$, ( $Q=\min (n_2-n,n-n_1),(varieswithn)$ $and$ $p\in [-Q,Q],n\in [n_1,n_2]$ )

令

$n{\mathrm{\Delta }}_t$

共有

$T$

點，

$m{\mathrm{\Delta }}_f$

共有

$F$

點，此算法其複雜度為

$\theta (TFmean(2Q+1))$

$(mean(Q)\cong \frac{n_2-n_1}{4},mean(2Q+1)\cong \frac{n_2-n_1}{2}+1\cong \frac{T}{2}(\because T=n_2-n_1+1))$

$\Rightarrow \theta (TFmean(2Q+1))=\theta (TF\frac{T}{2})=\theta (T^{2}F)$

離散傅立葉轉換(Using Discrete time Fourier Transform)

當 ${\mathrm{\Delta }}_t{\mathrm{\Delta }}_f=\frac{1}{2N},andN\ge 2Q+1$,

$W_x(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=2{\mathrm{\Delta }}_t\sum _{p=-Q}^{Q}x((n+p){\mathrm{\Delta }}_t)x^{*}((n-p){\mathrm{\Delta }}_t)e^{-j\frac{2\pi mp}{N}}$, 令 $q=p+Q\to p=q-Q$

$\Rightarrow W_x(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=2{\mathrm{\Delta }}_t{e}^{j\frac{2\pi mQ}{N}}\sum _{q=0}^{2Q}x((n+q-Q){\mathrm{\Delta }}_t)x^{*}((n-q+Q){\mathrm{\Delta }}_t)e^{-j\frac{2\pi mq}{N}}$

$\Rightarrow W_x(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=2{\mathrm{\Delta }}_t{e}^{j\frac{2\pi mQ}{N}}\sum _{q=0}^{N-1}{c}_1(q)e^{-j\frac{2\pi mq}{N}}$

其中：

$c_1(q)=x((n+q-Q){\mathrm{\Delta }}_t)x^{*}((n-q+Q){\mathrm{\Delta }}_t),for0\le q\le 2Q$

$c_1(q)=0,for2Q+1\le q\le N-1$

此實現方式的複雜度為 $\theta (TN\log (N))(\because FourierTransform:N-points\to N\log (N))$

Chirp-Z轉換

$W_x(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=2\sum _{p=-Q}^{Q}x((n+p){\mathrm{\Delta }}_t)x^{*}((n-p){\mathrm{\Delta }}_t)e^{-j4\pi mp{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{\mathrm{\Delta }}_t$

$\Rightarrow W_x(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=2{\mathrm{\Delta }}_t{e}^{-j2\pi m^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}\sum _{p=-Q}^{Q}x((n+p){\mathrm{\Delta }}_t)x^{*}((n-p){\mathrm{\Delta }}_t)e^{-j2\pi p^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{e}^{j2\pi (p-m{)}^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}$

Step1. $x_1(n,p)=x((n+p){\mathrm{\Delta }}_t)x^{*}((n-p){\mathrm{\Delta }}_t)e^{-j2\pi p^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}$

Step2. $X_2[n,m]=\sum _{p=-Q}^Q{x}_1[n,p]c[m-p],c[m]=e^{j2\pi m^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}$

Step3. $X(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=2{\mathrm{\Delta }}_t{e}^{-j2\pi m^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{X}_2[n,m]$

此實現方式複雜度一樣為 $\theta (TN\log (N))$ 但相較使用離散傅立葉轉換方式而言( $[\theta (T3N\log (N))=\theta (TN\log (N))]$ 數字3來自於convolution運算 IFT(FT * FT) )，速度來的更快。

資料來源

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