基于改进的Wigner-Ville分布与Chirp-Z变换的高分辨时频分析方法

High-Resolution Time-Frequency Analysis Based on Improved Wigner-Ville Distribution with Chirp-Z Transformation

摘要: 以Wigner-Ville分布(WVD)为代表的二次非线性时频分析方法具有良好的时频聚焦性，适合用于非平稳信号的时频分析，但在处理复杂多分量信号时，会产生严重的交叉干扰项影响人们对信号频谱的分析判断。伪Wigner-Ville分布(PWVD)和平滑伪Wigner-Ville分布(SPWVD)在WVD的基础上加窗改进核函数，可以有效抑制交叉项干扰，但也会降低对信号的时频分辨率。本文使用Chirp-Z变换(CZT)替换PWVD和SPWVD中的快速傅里叶变换(FFT)，得到两种高分辨时频分析方法PWVD-CZT和SPWVD-CZT，通过数值实验验证两种方法均可有效提高对信号的时频分辨率，SPWVD-CZT处理主频低、频带窄的复杂多分量非平稳信号具有更大的优势。

Abstract: The quadratic nonlinear time-frequency analysis method represented by the Wigner-Ville distribution (WVD) has good time-frequency focusing and is suitable for the time-frequency analysis of non-smooth signals, but when dealing with the complex multi-component signals, it will produce serious cross-interference terms affecting people’s analytical judgment of the signal spectrum. Pseudo Wigner-Ville distribution (PWVD) and smoothed pseudo Wigner-Ville distribution (SPWVD) improve the kernel function by adding a window on top of WVD, which can effectively suppress the cross term interference, but also reduce the time-frequency resolution of the signal. In this paper, the Chirp-Z transform (CZT) is used to replace the fast Fourier transform (FFT) in PWVD and SPWVD to obtain two high-resolution time-frequency analysis methods, PWVD-CZT and SPWVD-CZT, and it is verified through numerical experiments that the two methods are effective in improving the time-frequency resolution of the signal, and SPWVD-CZT has a greater advantage in dealing with the complex multi-component non-stationary signals with low primary frequency and narrow frequency band. The SPWVD-CZT has more advantages in dealing with complex multi-component non-stationary signals with low primary frequency and narrow frequency bands.

1. 引言

在信号处理领域，可以通过对信号应用时频分析尽可能多的提取信号中的信息，有助于人们解决实际工程类问题[1]。傅里叶变换可以将信号从时间域变换到频率域，从而对信号进行频域分析，包括幅值谱分析、相位谱分析和功率谱分析[2] [3]，与平稳信号不同，非平稳信号的时域或频域特性会随时间发生变化，傅里叶变换在整体上将信号分解为不同的频率分量，缺乏局域性信息，分析非平稳信号存在局限性[4]。

时频分析的研究始于20世纪40年代，1932年，Wigner首次提出了分辨率高，数学特性良好的Wigner分布，并将这一方法应用于量子力学领域[5]。1946年，Gabor将一维信号表示在时间和频率的二维坐标中，并推导出了时间-频率不确定性原理，将分析信号引入时频分析，提出了Gabor变换[6]，它可以反映信号在时域和频域的联合分布信息，为时频分析在信号分析中的应用奠定了基础。1948年，Ville将其应用于信号分析领域，推导出统计量子力学分布中的Wigner分布[7]，后发展为二次非线性时频Wigner-Ville分布，这种方法在处理单分量信号时具有时移不变性和良好的时频分辨率，但在处理多分量信号时，由于会产生交叉干扰项，致使时频分辨率较差。1966年，Cohen发现大多数时频分析方法可以通过对Wigner-Ville分布的核函数加权进行平滑处理得到，诞生了著名的Cohen类时频分布[8]，可以改善交叉项的干扰，从而提高信号分量的识别分析效果。

时频分析是一种在时域和频域联合分析和处理信号的方法，其中信号被表示为时间和频率的联合函数，并可以用信号的频谱成分作为时间函数的图像来表示。Wigner-Ville分布(WVD)是一种二次非线性时频分布，其在处理单分量信号时具有良好的时频聚焦性，但在处理多分量信号时会产生交叉项，与信号项不同，交叉项是信号时频分布里的干扰产物，它们往往展现出与原信号的物理性质相矛盾的结果[9]。早在1984年，有学者系统研究了WVD的交叉项问题，并提出了伪Wigner-Ville分布(PWVD)的概念，通过引入平滑窗函数抑制交叉项。为了进一步抑制交叉项，学者们在PWVD的基础上引入了频域平滑窗函数，提出了平滑伪Wigner-Ville分布(SPWVD)，SPWVD通过在时域和频域同时引入平滑窗函数，显著减少了交叉项干扰[10]-[12]。通过改进Wigner-Ville分布的核函数可以实现交叉项的抑制，伪Wigner-Ville分布(PWVD)和平滑伪Wigner-Ville分布(SPWVD)对交叉项具有良好的抑制效果，可以减少多分量信号时频分布的干扰因素，但同时也会降低对信号的时频聚焦性和时频分辨率。

随着研究的深入，有学者通过动态调整窗函数宽度，平衡时频分辨率和交叉项抑制[13]，基于信号瞬时频率提出了自适应窗函数设计方法，显著提高了PWVD的时频分辨率[14]，这些研究虽然可以在一定程度上提升时频分辨率，但过于依赖窗函数的选择，缺乏通用性，在处理多分量信号时，交叉项抑制效果有限。针对多分量信号处理方法，有学者研究了窗函数优化问题，提出了基于信号特性的自适应平滑方法[15]，但对于高度非平稳的信号，自适应平滑方法可能无法有效跟踪信号的时频变化。

针对复杂多分量非平稳信号，本文提出以Chirp-Z变换(CZT)替换伪Wigner-Ville分布和平滑伪Wigner-Ville分布中的快速傅里叶变换(FFT)，得到PWVD-CZT和SPWVD-CZT高分辨时频分析方法，通过调整采样间隔和采样点数实现对数据的插值计算，在原有方法的基础上进一步细化实频谱图，提高对信号的时频分辨率，为小目标信号的分析识别提供更加精确的方法支撑。我们分别使用生成的雷克子波和实际地震波信号进行数值实验，结果表明改进后的时频分析方法能够有效抑制多分量信号交叉项干扰，且对信号频谱局部细化，提高时频分辨率，是人们分析处理复杂多分量非平稳信号的有力工具。

本文结构安排如下：第一部分为理论部分，介绍改进的Wigner-Ville分布和Chirp-Z变换的定义和特性，阐述本文提出的时频分析方法PWVD-CZT和SPWVD-CZT的实现过程；第二部分是数值实验部分，使用地震信号验证本文所提方法的时频分析效果，并对比PWVD-CZT和SPWVD-CZT的效果差异；最后是结论部分。

2. 方法原理

2.1. PWVD和SPWVD

Wigner-Ville分布是最基本的二次非线性时频分析方法，也是在处理非平稳信号领域应用最多的时频分布方法。定义信号$x\left(t\right)$ 的Wigner-Ville分布如下：

${\text{WVD}}_x\left(t,f\right)={\displaystyle {\int }_{-\infty }^{+\infty }x\left(t+\frac{\tau }{2}\right)x^{*}\left(t-\frac{\tau }{2}\right){\text{e}}^{-j2\pi f\tau }\text{d}\tau }$ (1)

其中$t$ 是时间变量，$f$ 是频率变量，$\tau$ 是时间延迟，$x^{*}$ 是$x$ 的共轭，$j$ 是虚数单位(下同)。WVD未施加任何窗函数，而是将解析信号$x\left(t\right)$ 的逆向时移作为窗口，具有双线性的特征。

设观测信号包括两个源信号成分，$x=s_1+s_2$ ，那么：

$\begin{array}{c}{W}_x\left(t,f\right)={\displaystyle {\int }_{-\infty }^{+\infty }\left[s_1\left(t+\frac{\tau }{2}\right)+s_2\left(t+\frac{\tau }{2}\right)\right]}\left[s_1^{\ast }\left(t-\frac{\tau }{2}\right)+s_2^{\ast }\left(t-\frac{\tau }{2}\right)\right]{\text{e}}^{-j2\pi f\tau }\text{d}\tau \\ ={\displaystyle {\int }_{-\infty }^{+\infty }{s}_1}\left(t+\frac{\tau }{2}\right)s_1^{\ast }\left(t-\frac{\tau }{2}\right){\text{e}}^{-j2\pi f\tau }\text{d}\tau +{\displaystyle {\int }_{-\infty }^{+\infty }{s}_2}\left(t+\frac{\tau }{2}\right)s_2^{\ast }\left(t-\frac{\tau }{2}\right){\text{e}}^{-j2\pi f\tau }\text{d}\tau \\ +{\displaystyle {\int }_{-\infty }^{+\infty }\left[s_1\left(t+\frac{\tau }{2}\right)s_2^{\ast }\left(t-\frac{\tau }{2}\right)+s_2\left(t+\frac{\tau }{2}\right)s_1{}^{*}\left(t-\frac{\tau }{2}\right)\right]{\text{e}}^{-j2\pi f\tau }\text{d}\tau }\end{array}$ (2)

即：

$\text{WVD}\left(x;t,f\right)=\text{WVD}\left(s_1;t,f\right)+\text{WVD}\left(s_2;t,f\right)+\text{WVD}\left(s_1+s_2\right)$ (3)

WVD具有非线性的特征，多分量信号之和的WVD并不等于各分量信号的WVD之和，每两个信号项之间均会产生一个交叉项，干扰我们对各信号分量的正确判断[16]。因此，抑制交叉项成为二次型时频分布的重点研究课题。

图1是具有四个高斯分量信号的WVD时频分布图，可以看到每两个信号的能量响应之间有一些其他的能量分布，且会影响到主信号项成分的能量分布，这就是交叉干扰项。

Figure 1. WVD distribution of the four Gaussian component signals

图1. 四个高斯分量信号的WVD分布

由两分量信号推广到多分量信号的情况，则在信号的WVD时频分布中，共同组成信号表示的有各个主信号项分量，与此同时每两个主信号项分量二者之间还有一个交叉项分量存在。对于一个含有$p$ 个分量的信号$x\left(t\right)$ ，其${\text{WVD}}_x\left(t,f\right)$ 中将含有$p$ 个信号项成分和$\left(\begin{array}{c}p\\ 2\end{array}\right)=\frac{p\left(p-1\right)}{2}$ 个交叉项成分，增加了对多分量信号的时频分析难度。

通过改进WVD的核函数可以实现对交叉项的抑制，其中最简单的一种做法是通过对变量$\tau$ 施加窗函数$h\left(\tau \right)$ ，这相当于对原始信号在时间方向上施加了一个可滑动的时间窗口函数，可以抑制时间方向上的交叉项干扰，这样改进后的WVD称为伪Wigner-Ville分布(PWVD)：

${\text{PWVD}}_x\left(t,f\right)={\displaystyle {\int }_{-\infty }^{+\infty }x\left(t+\frac{\tau }{2}\right)x^{*}}\left(t-\frac{\tau }{2}\right)h\left(\tau \right){\text{e}}^{-j2\pi f\tau }\text{d}\tau$ (4)

图2是包含四个高斯分量信号的PWVD时频分布图，与图1对比发现，水平方向相邻两个主信号项之间的交叉干扰项被抑制了，这和PWVD分布在时间方向上施加了窗函数$h\left(\tau \right)$ 有关。

为进一步抑制频率方向上的交叉项干扰，继续对变量$u$ 施加一个窗函数$g\left(u\right)$ ，这样得到的改进的Wigner-Ville分布在时域和频域兼具卷积平滑作用，具备更加理想的抑制交叉项效果，称为平滑伪Wigner-Ville分布(SPWVD)：

Figure 2. PWVD distribution of the four Gaussian component signals

图2. 四个高斯分量信号的PWVD分布

${\text{SPWVD}}_x\left(t,f\right)={\displaystyle {\int }_{-\infty }^{+\infty }{\displaystyle {\int }_{-\infty }^{+\infty }g\left(u\right)h\left(\tau \right)x\left(t-u+\frac{\tau }{2}\right)x^{*}}}\left(t-u-\frac{\tau }{2}\right){\text{e}}^{-j2\pi f\tau }\text{d}u\text{d}\tau$ (5)

式中$h\left(\tau \right)$ 和$g\left(u\right)$ 是两个实的偶窗函数。

图3是包含同样四个高斯分量信号的SPWVD时频分布图，和图2对比可以看到垂直方向的每两个信号之间的交叉干扰项也得到了抑制，这是因为SPWVD在时间和频率方向上都施加了窗函数，对时间频率两方向都有拉伸平滑作用，具有较佳的交叉项抑制效果，同时保持较好的时频特性和聚焦性能。

2.2. Chirp-Z变换

离散傅里叶变换(DFT)本质上是在单位圆上对一个有限长序列的Z变换做等间隔采样，在实际应用中，DFT具有很大的局限性，一个有限长序列的N个点对应于均匀分布在Z平面单位圆周上的N个采样点，它的频率分辨率完全取决于信号序列的长度N，信号的频率分辨率受限于信号序列的长度。快速傅里叶变换(FFT)可以很快的计算出单位圆上等距的全部采样值，但有时并不是所有的采样点都具有意义[17]。Chirp-Z变换(CZT)可以沿着Z平面上不限于单位圆的路径并快速实现Z变换，突破了FFT算法对信号序列长度和路径的限制[18]。

设$X\left(z\right)$ 是一个长度为N的离散序列$x\left(n\right)$ 的Z变换：

$X\left(z\right)={\displaystyle \sum _{n=0}^{N-1}x\left(n\right)z^{-n}}$ (6)

$n$ 为离散时间变量，$z$ 为Z平面上的采样因子。

Figure 3. SPWVD distribution of the four Gaussian component signals

图3. 四个高斯分量信号的SPWVD分布

Z平面上$z_k$ 点的采样值为：

$z_k=A{W}^{-k},k=0,1,\cdots ,M-1$ (7)

$M$ 为Z变换的长度，$A$ 为起始位置。

起始位置可以用它的半径$A_0$ 和幅角${\theta }_0$ 表示，$A_0$ 可取任意正实数：

$A=A_0{\text{e}}^{j{\theta }_0}$ (8)

参数$W$ 可表示为：

$W=W_0{\text{e}}^{j{\varphi }_0}$ (9)

$W_0$ 为螺线的伸展率，可取任意正实数，当$W_0>1$ 时，随着$k$ 的增大螺线趋向圆心，当$W_0<1$ 时，随着$k$ 的增大螺线趋向圆外，在对信号作频谱分析时需要在单位圆上计算CZT，此时$W_0=1$ 。${\varphi }_0$ 是螺线上采样点之间的角度，由于${\varphi }_0$ 可以取任意的值，所以CZT算法的分辨率可以是任意的值，不受信号序列长度N的约束，应用起来更加灵活方便，有利于研究具有任意起始频率的高分辨率窄带频谱。

因此有：

$z_k=A_0{\text{e}}^{j{\theta }_0}{W}_0{\text{e}}^{-j{\varphi }_{0}k},k=0,1,\cdots ,M-1$ (10)

沿Z平面上采样点$z_k$ 的Z变换为：

$A\left(z_k\right)={\displaystyle \sum _{n=0}^{N-1}x\left(n\right)z^{-n}}={\displaystyle \sum _{n=0}^{N-1}x\left(n\right)A^{-n}{W}^{nk}},k=0,1,\cdots ,M-1$ (11)

为了提高运算速度，把乘积项换成求和项，由Bluestein公式：

$nk=\frac{n^2+k^2-{\left(n-k\right)}^2}{2},k=0,1,\cdots ,M-1$ (12)

于是，$z_k$ 的Z变换又可改写为：

$X\left(z_k\right)=W^{\frac{k^2}{2}}{\displaystyle \sum _{n=0}^{N-1}x\left(n\right)A^{-n}{W}^{\frac{n^2}{2}}{W}^{\frac{-{\left(n-k\right)}^2}{2}}},k=0,1,\cdots ,M-1$ (13)

令：

$q\left(n\right)=x\left(n\right)A^{-n}{W}^{\frac{n^2}{2}}$ (14)

$g\left(n\right)=W^{\frac{-n^2}{2}}$ (15)

则上式可简写为：

$X\left(z_k\right)=W^{\frac{k^2}{2}}{\displaystyle \sum _{n=0}^{N-1}q\left(n\right)g\left(n-k\right),k=0,1,\cdots ,M-1}$ (16)

考虑信号在Z平面上采样点$z_k$ 都在单位圆上的情况，这意味着Chirp-Z变换类似于傅里叶变换，不同于DFT需要对Z平面上的N点数据进行N点变换，CZT不必变换全部的点，只需考虑变换局部的点，因此不需要使$M=N$ 。

2.3. PWVD-CZT和SPWVD-CZT

PWVD在WVD的基础上，通过给核函数加权，对信号在时间方向上施加一个可移动的窗函数，可以在很大程度上抑制信号项成分之间交叉项干扰的产生。CZT算法通过调整采样间隔和采样点数，可以实现时频谱图的细化，提升信号的时频分辨率。将PWVD和CZT结合，得到高分辨时频分析方法PWVD-CZT，实现过程如下：

首先对长度为$N$ 的离散信号序列$x\left(n\right)$ 进行Hilbert变换得到解析信号$z\left(n\right)$ ，取解析信号的复模量会产生包络效应，可以消除波形中的震荡。

离散信号序列$x\left(n\right)$ 的解析信号$z\left(n\right)$ 为：

$z\left(n\right)=x\left(n\right)+jH\left[x\left(n\right)\right]$ (17)

$x\left(n\right)$ 是解析信号$z\left(n\right)$ 的实部，$x\left(n\right)$ 的Hilbert变换$H\left[x\left(n\right)\right]$ 是$z\left(n\right)$ 的虚部。

解析信号$z\left(n\right)$ 的WVD分布为：

$\text{WVD}\left(n,\sigma \right)=2{\displaystyle \sum _{-\infty }^{+\infty }z\left(n+l\right)z^{*}\left(n-l\right){\text{e}}^{-j2l\sigma }}$ (18)

$\text{WVD}\left(n,\sigma \right)$ 是解析信号$z\left(n\right)$ 的Wigner-Ville分布，$l$ 为时间步长，$z^{*}$ 是$z$ 的共轭，${\text{e}}^{-j2l\sigma }$ 是傅里叶变换因子，$\sigma$ 是频率因子。

解析信号$z\left(n\right)$ 的瞬时自相关函数为：

$r\left(n,l\right)=z\left(n,l\right)z{\left(n-l\right)}^{*}$ (19)

对瞬时自相关函数$r\left(n,l\right)$ 乘一个长度为$\left(\lfloor \frac{N}{4}\rfloor +1\right)$ 的汉明窗$h\left(l\right)$ ，且当$\left|l\right|>\lfloor \frac{N}{8}\rfloor$ 时，$h\left(l\right)=0$ ，得到加窗后的离散瞬时自相关函数为：

$R\left(n,l\right)=z\left(n+l\right)z{\left(n-l\right)}^{*}h\left(l\right)h{\left(-l\right)}^{*}$ (20)

$h{\left(-l\right)}^{*}$ 是窗函数$h\left(l\right)$ 的共轭。

令：

$Z\left(n,l\right)=z\left(n+l\right)h\left(l\right)$ (21)

得到PWVD的核函数为：

$K_n\left(l\right)=\{\begin{array}{ll}Z\left(n,l\right)Z{\left(n,-l\right)}^{*},\hfill & \left|l\right|\le \min \left\{n,N-n\right\}\hfill \\ 0,\hfill & \left|l\right|>\min \left\{n,N-n\right\}\hfill \end{array}$ (22)

对核函数$K_n\left(l\right)$ 做循环位移，使得$l$ 的范围位于$0~\left(N-1\right)$ 内，得到循环位移后的核函数$K_N\left(l\right)$ 为：

$K_N\left(l\right)=\{\begin{array}{ll}Z\left(n,l\right)Z{\left(n,-l\right)}^{*},\hfill & 0\le l\le \frac{N}{2}-1\hfill \\ Z\left(n,l-N\right)Z{\left(n,-l+N\right)}^{*},\hfill & \frac{N}{2}\le l\le N\hfill \end{array}$ (23)

$Z{\left(n,-l+N\right)}^{*}$ 是$Z\left(n,-l+N\right)$ 的共轭。

最后，用CZT替换PWVD算法中的FFT。对PWVD的核函数$K_N\left(l\right)$ 乘上加权系数$A^{-l}{W}^{\frac{l^2}{2}}$ ，得到序列$f^{*}\left(l\right)$ ：

$f^{*}\left(l\right)=K_N\left(l\right)A^{-l}{W}^{\frac{l^2}{2}},0\le l\le N-1$ (24)

由于WVD的周期为$\pi$ ，为了使得信号在Z平面上沿单位圆的采样点在整个单位圆周上均匀分布，因此CZT中的相角${\theta }_0$ 和幅角${\phi }_0$ 应扩大为原来的两倍，故有：

$A=A_0{\text{e}}^{j2{\theta }_0}$ (25)

$W=W_0{\text{e}}^{j2{\phi }_0}$ (26)

对序列$f^{*}\left(l\right)$ 做循环位移，得到序列$f\left(l\right)$ ：

$f\left(l\right)=\{\begin{array}{ll}{f}^{*}\left(l\right),\hfill & 0\le l\le \frac{N}{2}-1\hfill \\ 0,\hfill & \frac{N}{2}\le l\le Q-\frac{N}{2}\hfill \\ f^{*}\left(l-Q+N\right),\hfill & Q-\frac{N}{2}+1\le l\le Q-1\hfill \end{array}$ (27)

构造单位相应序列$\delta \left(l\right)$ ：

$\delta \left(l\right)=\{\begin{array}{ll}{W}^{\frac{l^2}{2}},\hfill & 0\le l\le M-1\hfill \\ 0,\hfill & M\le l\le Q-N\hfill \\ W^{-\frac{{\left(\frac{N}{2}-l\right)}^2}{2}},\hfill & Q-N+1\le l\le Q-1\hfill \end{array}$ (28)

$Q$ 是为求CZT构造的序列长度，其取值为最接近$\left(N+M-1\right)$ 的二次幂；$M$ 是CZT上的采样点数。

对$f\left(l\right)$ 和$\delta \left(l\right)$ 做傅里叶变换：

$F\left(r\right)=\text{FFT}\left[f\left(l\right)\right]$ (29)

$\Delta \left(r\right)=\text{FFT}\left[\delta \left(l\right)\right]$ (30)

将得到的$F\left(r\right)$ 和$\Delta \left(r\right)$ 在频域内相乘得到：

$Y\left(r\right)=F\left(r\right)\Delta \left(r\right)$ (31)

对$Y\left(r\right)$ 做傅里叶逆变换，得到：

$y\left(r\right)=\text{IFFT}\left[Y\left(r\right)\right]$ (32)

最后，对$y\left(r\right)$ 乘上加权系数，得到信号的PWVD-CZT分布结果：

$P_{\text{wvd-czt}}\left(z_l\right)=y\left(l\right)W^{\frac{l^2}{2}},0\le l\le M$ (33)

SPWVD-CZT算法的实现过程和PWVD-CZT算法的实现过程类似，区别在于对WVD的核函数进行加权的过程，PWVD-CZT算法仅在信号的时间方向乘上一个窗函数，SPWVD-CZT算法在信号的时间和频率方向分别乘上了一个窗函数，对信号在时频域均有平滑卷积效果，对多分量信号的交叉项抑制效果更优，但同时也增加了更多的计算空间复杂度。

对式(19)在频域方向上乘一个长度为$\left(\lfloor \frac{N}{4}\rfloor +1\right)$ 的汉明窗$g\left(l\right)$ ，且当$\left|l\right|>\lfloor \frac{N}{8}\rfloor$ 时，$g\left(l\right)=0$ ，得到加窗后的离散瞬时自相关函数为：

$R\left(n,l\right)=z\left(n+l\right)z{\left(n-l\right)}^{*}g\left(l\right)g{\left(-l\right)}^{*}$ (34)

$g{\left(-l\right)}^{*}$ 是窗函数$g\left(l\right)$ 的共轭。

再对瞬时自相关函数$r\left(n,l\right)$ 在时域方向上乘一个长度为$\left(\lfloor \frac{N}{10}\rfloor +1\right)$ 的汉明窗$h\left(l\right)$ ，且当$\left|l\right|>\lfloor \frac{N}{20}\rfloor$ 时，$h\left(l\right)=0$ ，得到加窗后的离散瞬时自相关函数为：

$R\left(n,l\right)=h\left(n-2l\right)z\left(n+l\right)z{\left(n-l\right)}^{*}g\left(l\right)g{\left(-l\right)}^{*}$ (35)

令：

$Z\left(n,l\right)=z\left(n+l\right)h\left(l\right)$ (36)

得到SPWVD的核函数为：

$K_n\left(l\right)=\{\begin{array}{ll}Z\left(n,l\right)Z{\left(n,-l\right)}^{*},\hfill & \left|l\right|\le \min \left\{n,N-n\right\}\hfill \\ 0,\hfill & \left|l\right|>\min \left\{n,N-n\right\}\hfill \end{array}$ (37)

经过式(23)~(32)的运算过程，最终得到信号的SPWVD-CZT分布结果：

$S_{\text{PWVD-CZT}}\left(z_l\right)=y\left(l\right)W^{\frac{l^2}{2}},0\le l\le M$ (38)

3. 数值实验

图4为来自某工区抽取的一道实际地震信号波形记录，采样频率为500 Hz，采样点数$N$ 为3002。本文使用这段地震信号数据分别对所提方法进行数值实验，对比实验结果。

Figure 4. Seismic signal waveform

图4. 地震信号波形图

图5(a)是图4展示的地震信号的PWVD频谱，图5(b)是地震信号的PWVD-CZT频谱，为清晰直观地对比两种方法对信号处理的效果，我们选取0.4~0.8 s频率为50~54 Hz的频段进行局部放大显示，图5(c)和图5(d)分别是图5(a)和图5(b)的局部放大图，从图中对比可以看出，PWVD-CZT方法比PWVD方法更能清晰地刻画出信号的能量分布与时间和频率的关系。

图6(a)是图4展示的地震信号的SPWVD频谱，图6(b)是地震信号的SPWVD-CZT频谱，为清晰直观地对比两种方法对信号处理的效果，我们选取0.4~1 s频率为54~58 Hz的频段进行局部放大显示，

Figure 5. (a) PWVD spectrum; (b) PWVD-CZT spectrum; (c) Fig. (a) local zoomed-in display; (d) Fig. (b) local zoomed-in display

图5. (a) PWVD频谱；(b) PWVD-CZT频谱；(c)图(a)局部放大显示；(d)图(b)局部放大显示

图6(c)和图6(d)分别是图6(a)和图6(b)的局部放大图，对比可以看出，地震信号的SPWVD-CZT时频分布细化了局部细节，具有更高的时频分辨率。图5(a)和图6(a)对比可以看出，在处理复杂地震信号时，对比PWVD，SPWVD可以保持更好的时频特性和聚焦性能，在频带较窄的情况下，信号的SPWVD频谱具有更高的时频分辨率。

Figure 6. (a) SPWVD spectrum; (b) SPWVD-CZT spectrum; (c) Fig. (a) local zoomed-in display; (d) Fig.; (b) local zoomed-in display

图6. (a) SPWVD频谱；(b) SPWVD-CZT频谱；(c)图(a)局部放大显示；(d)图；(b)局部放大显示

4. 结论

Wigner-Ville分布对信号具有良好的时频聚焦性，但在处理复杂多分量非平稳信号时会产生明显的交叉项干扰，通过施加窗函数改进核函数可以有效抑制交叉项，但也会降低对信号的时频分辨率。本文基于加窗后的Wigner-Ville分布，用CZT替换其中的FFT，得到PWVD-CZT和SPWVD-CZT两种时频分析方法，通过数值实验，验证了这两种时频分析方法可以有效提高信号的时频分辨率，其中SPWVD-CZT处理主频低、频带窄的复杂多分量非平稳信号具有更强的优势。

NOTES

^*通讯作者。