韋格納分佈的時頻分析

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

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

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

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

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

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

於是，通常我們會假設 $x(t)=0\ for\ t<n_{1}\Delta _{t}\ and\ t>n_{2}\Delta _{t}$，也就是 $x(t)$ 為一個有限長度的訊號：

如圖所示：

此時 $x((n+p)\Delta _{t})x^{*}((n-p)\Delta _{t})=0,\ if\ (n+p)\not \in [n_{1},n_{2}]\ or\ (n-p)\not \in [n1,n2]$

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

1. 對於$(n+p)$：$n_{1}\leq (n+p)\leq {n_{2}}\longrightarrow (n_{1}-n)\leq p\leq (n_{2}-n)$
2. 對於$(n-p)$：$n_{1}\leq (n-p)\leq {n_{2}}\longrightarrow (n_{1}-n)\leq -p\leq (n_{2}-n)$，也就是 $(n-n_{2})\leq p\leq (n-n_{1})$

於是便得知 $p$ 的範圍為：$\max {(n_{1}-n,n-n_{2})}\leq p\leq \min(n_{2}-n,n-n_{1})\longrightarrow -\min {(n_{2}-n,n-n_{1})}\leq p\leq \min(n_{2}-n,n-n_{1})\longrightarrow -Q\leq p\leq Q,\ Q=min(n_{2}-n,n-n_{1})$

用圖片來理解：

注意：當 $n>n_{2}$ 或 $n<n_{1}$時，找不到適當的 $p$ 來滿足不等式。

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

$if\ x(t)=0\ for\ t<n_{1}\Delta _{t}\ and\ t>n_{2}\Delta _{t}$，

$W_{x}(n\Delta _{t},m\Delta _{f})=2\sum _{p=-Q}^{Q}x((n+p)\Delta _{t})x^{*}((n-p)\Delta _{t})e^{-j4\pi mp\Delta _{t}\Delta _{f}}\Delta _{t}$

( $Q=\min {(n_{2}-n,n-n_{1})},(varies\ with\ n)$ $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\Delta _{t},m\Delta _{f})=2\sum _{p=-Q}^{Q}x((n+p)\Delta _{t})x^{*}((n-p)\Delta _{t})e^{-j4\pi mp\Delta _{t}\Delta _{f}}\Delta _{t}$, ( $Q=\min {(n_{2}-n,n-n_{1})},(varies\ with\ n)$ $and$ $p\in [-Q,Q],\ n\in [n_{1},n_{2}]$ )

令$n\Delta _{t}$共有$T$點，$m\Delta _{f}$共有$F$點，此算法其複雜度為 $\theta (TF\ mean(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 (TF\ mean(2Q+1))=\theta (TF{\frac {T}{2}})=\theta {(T^{2}F)}$

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

當 $\Delta _{t}\Delta _{f}={\frac {1}{2N}},\ and\ N\geq 2Q+1$,

$W_{x}(n\Delta _{t},m\Delta _{f})=2\Delta _{t}\sum _{p=-Q}^{Q}x((n+p)\Delta _{t})x^{*}((n-p)\Delta _{t})e^{-j{\frac {2\pi mp}{N}}}$, 令 $q=p+Q\rightarrow p=q-Q$

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

$\Rightarrow W_{x}(n\Delta _{t},m\Delta _{f})=2\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)\Delta _{t})x^{*}((n-q+Q)\Delta _{t}),\ for\ 0\leq q\leq 2Q$

$c_{1}(q)=0,\ for\ 2Q+1\leq q\leq N-1$

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

Chirp-Z轉換

$W_{x}(n\Delta _{t},m\Delta _{f})=2\sum _{p=-Q}^{Q}x((n+p)\Delta _{t})x^{*}((n-p)\Delta _{t})e^{-j4\pi mp\Delta _{t}\Delta _{f}}\Delta _{t}$

$\Rightarrow W_{x}(n\Delta _{t},m\Delta _{f})=2\Delta _{t}e^{-j2\pi m^{2}\Delta _{t}\Delta _{f}}\sum _{p=-Q}^{Q}x((n+p)\Delta _{t})x^{*}((n-p)\Delta _{t})e^{-j2\pi p^{2}\Delta _{t}\Delta _{f}}e^{j2\pi (p-m)^{2}\Delta _{t}\Delta _{f}}$

Step1. $x_{1}(n,p)=x((n+p)\Delta _{t})x^{*}((n-p)\Delta _{t})e^{-j2\pi p^{2}\Delta _{t}\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}\Delta _{t}\Delta _{f}}$

Step3. $X(n\Delta _{t},m\Delta _{f})=2\Delta _{t}e^{-j2\pi m^{2}\Delta _{t}\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.