韦格纳分布的时频分析

韦格纳分布Wigner Distribution Function(缩写为WDF)，是时频分析中的一种分析方式，以下为其方程式的转换：

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

经整理后得到 $W_x(t,f)=2{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}x(t+{\tau }^{\prime })\cdot x^{\ast }(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^{\ast }((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^{\ast }((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^{\ast }((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^{\ast }((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^{\ast }((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^{\ast }((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^{\ast }((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^{\ast }((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^{\ast }((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^{\ast }((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.