查看“︁韋格納分佈的時頻分析”︁的源代码



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

<math>W_x(t, f) = \int_{-\infty}^{\infty} x(t + \tau/2) \cdot{x^*(t-\tau/2)e^{-j2\pi f}\cdot{d\tau}}</math>

經整理後得到  <math>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</math>

連續訊號要得到其數位實現，必須經過採樣(sampling), 使得 <math>t = n \Delta_t, \ f = m \Delta_f, \tau'= p \Delta_t</math>

重新回顧上述式子可以得到，<math>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}</math>

當 <math>x(t)</math> 不是 time-limited signal，會很難實現。


於是，通常我們會假設 <math>x(t) = 0 \ for \ t < n_1 \Delta_t \  and \ t > n_2 \Delta_t</math>，也就是 <math>x(t)</math> 為一個有限長度的訊號：

如圖所示：

[[File:X(t)_discrete_signal.png|替代=|无框|459x459像素]]

此時 <math>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]</math>

繼續討論 <math>p</math> 的範圍 ( 於<math>n</math>為一固定值時 )：

# 對於<math>(n + p)</math>：<math>n_1 \leq (n+p) \leq {n_2}\longrightarrow (n_1-n) \leq p \leq (n_2-n)</math>
# 對於<math>(n - p)</math>：<math>n_1 \leq (n-p) \leq {n_2}\longrightarrow (n_1-n) \leq -p \leq (n_2-n)</math>，也就是 <math>(n-n_2) \leq p \leq (n-n_1)</math>

於是便得知 <math>p</math> 的範圍為：<math>\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)</math>

用圖片來理解：

[[File:X(t)_delta.png|替代=|无框|455x455像素]]

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

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

<math>if \ x(t) = 0 \ for \ t < n_1\Delta_t \ and \ t > n_2\Delta_t</math>，

<math>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</math><blockquote>(   <math>Q = \min{(n_2 - n, n - n_1)}, (varies \ with \ n)</math> <math>and</math>  <math>p \in [-Q, Q], \ n \in [n_1, n_2]</math>   )</blockquote>以下提供3種實現方式：

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

<br />
==== 暴力法(Direct Implementation) ====
根據 <math>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</math>, (   <math>Q = \min{(n_2 - n, n - n_1)}, (varies \ with \ n)</math> <math>and</math>  <math>p \in [-Q, Q], \ n \in [n_1, n_2]</math>   )

令<math>n\Delta_t</math>共有<math>T</math>點，<math>m\Delta_f</math>共有<math>F</math>點，此算法其複雜度為 <math>\theta(TF \ mean(2Q+1))</math><blockquote>[[File:WDF_Q_Direct_implement.png|替代=|无框]]

<math>( \ 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) \ )</math></blockquote><math>\Rightarrow \theta(TF \ mean(2Q+1)) = \theta(TF\frac{T}{2}) = \theta{(T^2F)}</math>
<br />

==== 離散傅立葉轉換(Using Discrete time Fourier Transform) ====
當 <math>\Delta_t \Delta_f = \frac{1}{2N}, \ and \ N \geq 2Q + 1</math>, 

<math>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}}</math>, 令 <math>q = p + Q \rightarrow p = q - Q</math>

<math>\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}} </math>

<math>\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}}
 </math>

其中：

<math>c_1(q) = x((n+q-Q)\Delta_t)x^*((n-q+Q)\Delta_t), \ for\ 0 \leq q \leq 2Q
 </math>

<math>c_1(q) = 0, \ for\ 2Q+1 \leq q \leq N-1
 </math>

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

<br />

==== Chirp-Z轉換 ====
<math>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</math>

<math>\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}</math>

Step1. <math>x_1(n, p) = x((n+p)\Delta_t)x^*((n-p)\Delta_t)e^{-j2\pi p^2\Delta_t\Delta_f}</math>

Step2. <math>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}</math>

Step3. <math>X(n\Delta_t, m\Delta_f) = 2\Delta_t e^{-j2\pi m^2\Delta_t\Delta_f}X_2[n, m]</math>


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

<br />

=== 資料來源 ===

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

<br />

[[Category:時頻分析]]