论文阅读——复频谱的局部化：S变换

复频谱的局部化：S变换

R. G. Stockwell, L. Mansinha and R. P. Lowe, “Localization of the complex spectrum: the S transform,” in IEEE Transactions on Signal Processing, vol. 44, no. 4, pp. 998-1001, April 1996, doi: 10.1109/78.492555.

引言

在地球物理数据分析以及许多其他学科中，平稳时间序列的概念是一种数学理想化，它从未真正实现，在信号到达检测中也不是特别有用。虽然整个时间序列的傅里叶变换确实包含了时间序列中频谱成分的信息，但对于大量实际应用而言，这些信息是不充分的。

以地震学中的地震记录为例来说明这个问题的重要性：从地震传来的第一个信号是P波（初波），随后是沿不同路径传播的其他P波。P波到达后紧接着是S波（次波）和更高振幅的频散表面波。这些振荡的振幅可以在P波到达后的几分钟内增加两个数量级以上。这种时间序列的频谱成分显然具有很强的时间依赖性。如果仅使用传统的傅里叶变换，我们只能得到整个信号的频率成分，而无法知道这些频率成分在时间上的分布。因此，拥有一种联合时频表示（TFR）是非常必要的。

本文提出了一种新的变换——S变换，它提供了一种具有频率依赖分辨率的时频表示，同时通过时间平均与傅里叶谱保持直接关系。过去已经提出了几种检查频谱时变特性的技术，其中包括Gabor变换、相关的短时傅里叶变换（STFT）、连续小波变换（CWT）以及称为Cohen类的双线性时频分布类，其中Wigner分布是其成员之一。

S变换的理论基础

从连续小波变换到S变换

S变换可以通过多种方法得出。我们认为将S变换推导为连续小波变换（CWT）的"相位校正"是很有启发性的。函数$h(t)$的连续小波变换$W(\tau ,d)$定义为：

$$W(\tau ,d)={\int }_{-\infty }^{\infty }h(t)w(t-\tau ,d)dt$$

其中$w(t,d)$是基本母小波的缩放副本，膨胀因子$d$决定了小波$w(t,d)$的"宽度"，从而控制分辨率。与此同时，存在一个关于母小波$w(t,d)$的可容许条件，即$w(t,d)$必须具有零均值。

函数$h(t)$的S变换定义为具有特定母小波的CWT乘以相位因子：

$$S(\tau ,f)=e^{i2\pi f\tau }W(\tau ,d)$$

其中母小波定义为：

$$w(t,f)=\frac{|f|}{\sqrt{2\pi }}{e}^{-\frac{t^2{f}^2}{2}}{e}^{-i2\pi ft}$$

这里需要注意的是，膨胀因子$d$是频率$f$的倒数。方程(3)中的小波不满足可容许小波的零均值条件，因此(2)严格来说不是CWT。明确写出，S变换为：

$$S(\tau ,f)={\int }_{-\infty }^{\infty }h(t)\frac{|f|}{\sqrt{2\pi }}{e}^{-\frac{(t-\tau {)}^2{f}^2}{2}}{e}^{i2\pi ft}dt$$

这个表达式清楚地显示了S变换的核心思想：使用一个高斯窗口来局部化信号，窗口的宽度与频率成反比（频率越高，窗口越窄），同时保持调制正弦函数相对于时间轴固定。

S变换与傅里叶变换的关系

如果S变换确实是局部频谱的表示，人们期望通过简单的时间平均操作就能得到傅里叶谱。事实上，容易证明：

$${\int }_{-\infty }^{\infty }S(\tau ,f)d\tau =H(f)$$

其中$H(f)$是$h(t)$的傅里叶变换。这个性质表明$h(t)$可以从$S(\tau ,f)$精确恢复：

$$h(t)={\int }_{-\infty }^{\infty }\left\{{\int }_{-\infty }^{\infty }S(\tau ,f)d\tau \right\}{e}^{i2\pi ft}df$$

这个恢复公式显然不同于CWT的概念，它确保了S变换的可逆性。

瞬时频率的推广

S变换提供了瞬时频率（IF）概念向宽带信号的扩展。变量$\tau$和固定参数$f_1$定义的一维函数$S(\tau ,f_1)$称为一个声音（voice，这是从CWT借用的术语）。特定频率$f_1$的声音可以写成：

$$S(\tau ,f_1)=A(\tau ,f_1)e^{i\Phi (\tau ,f_1)}$$

由于声音隔离了特定成分，可以使用方程(7)中的相位来确定Bracewell定义的瞬时频率：

IF(τ,f1)=12π∂∂τ{2πf1τ+Φ(τ,f1)}IF(\tau, f_1) = \frac{1}{2\pi}\frac{\partial}{\partial \tau}\{2\pi f_1\tau + \Phi(\tau, f_1)\}IF(τ,f1​)=2π1​∂τ∂​{2πf1​τ+Φ(τ,f1​)}

这个定义的有效性可以通过简单的例子来验证。对于$h(t)=\cos (2\pi \omega t)$的情况，函数$\Phi (\tau ,f)=2\pi (\omega -f)\tau$，瞬时频率恰好等于$\omega$。

线性性质与噪声处理

S变换的线性性质确保对于加性噪声的情况（可以将数据建模为data(t) = signal(t) + noise(t)），S变换给出：

S{data}=S{signal}+S{noise}S\{data\} = S\{signal\} + S\{noise\}S{data}=S{signal}+S{noise}

这是相对于Cohen类双线性时频分布的优势，在双线性类中：

TFR{data}=TFR{signal}+TFR{noise}+2×TFR{signal}×TFR{noise}TFR\{data\} = TFR\{signal\} + TFR\{noise\} + 2 \times TFR\{signal\} \times TFR\{noise\}TFR{data}=TFR{signal}+TFR{noise}+2×TFR{signal}×TFR{noise}

交叉项的存在使得在噪声局部谱上建立可靠的信号估计变得困难。

离散S变换的实现

离散化过程

设$h[kT]$，$k=0,1,...,N-1$表示对应于$h(t)$的离散时间序列，时间采样间隔为$T$。离散傅里叶变换由下式给出：

$$H\left[\frac{n}{NT}\right]=\frac{1}{N}\sum _{k=0}^{N-1}h[kT]e^{-\frac{i2\pi nk}{N}}$$

其中$n=0,1,...,N-1$。在离散情况下，S变换是时间序列定义的向量在一组张成向量上的投影。这些张成向量不是正交的，S变换的元素也不是独立的。每个基向量（傅里叶变换的）通过与$N$个移位高斯函数的逐元素乘积被分成$N$个局部化向量，使得这$N$个局部化向量的和是原始基向量。

使用方程(9)，离散时间序列$h[kT]$的S变换由下式给出（令$f\to n/NT$和$\tau \to jT$）：

$$S\left[jT,\frac{n}{NT}\right]=\sum _{m=0}^{N-1}H\left[\frac{m+n}{NT}\right]e^{-\frac{2{\pi }^2{m}^2}{n^2}}{e}^{\frac{i2\pi mj}{N}},n\ne 0$$

对于$n=0$的声音，它等于定义为以下的常数：

$$S[jT,0]=\frac{1}{N}\sum _{m=0}^{N-1}h\left(\frac{m}{NT}\right)$$

方程(12)将时间序列的常数平均值放入零频率声音中，从而确保逆变换是精确的。

离散逆S变换

离散S变换的逆变换由下式给出：

$$h[kT]=\sum _{n=0}^{N-1}\left\{\frac{1}{N}\sum _{j=0}^{N-1}S\left[jT,\frac{n}{NT}\right]\right\}{e}^{\frac{i2\pi nk}{N}}$$

实验结果与分析

图1：频率分辨率对比分析

图1(a) 显示了一个合成时间序列，由三个不同的成分组成：

• 前半部分（0-63采样点）：低频信号 $h[0:63]=\cos (2\pi t\times 6.0/128.0)$
• 后半部分（63-127采样点）：中频信号 $h[63:127]=\cos (2\pi t\times 25.0/128.0)$
• 高频突发（20-30采样点）：$h[20:30]=h[20:30]+0.5\cos (2\pi t\times 52.0/128.0)$

图1(b) 展示了S变换的振幅谱。可以清楚地看到S变换的频率依赖分辨率特性：

• 对于低频信号（$f=6/128$），S变换使用较宽的时间窗口，提供了良好的频率分辨率
• 对于高频突发（$f=52/128$），S变换使用较窄的时间窗口，精确地定位了突发的时间位置
• 中频信号在时频平面上也得到了清晰的表示

图1© 显示了使用固定高斯窗口（标准差=16单位）的短时傅里叶变换（STFT）。由于固定的窗口宽度，STFT无法同时提供良好的时间和频率分辨率。高频突发被时间窗口平滑，导致时间定位不准确。

图1(d) 展示了使用矩形窗（长度=20单位）的STFT结果。虽然使用了不同的窗口函数，但固定窗口宽度的根本限制仍然存在。

图2：交叉啁啾和突发信号检测

图2(a) 显示了一个更复杂的合成时间序列，包含：

• 两个交叉啁啾信号：
  • 上升啁啾：$\cos (2\pi (10+t/7)\times t/256)$
  • 下降啁啾：$\cos (2\pi (256/2.8-t/6.0)\times t/256)$
• 两个高频突发：
  • 第一个突发（114-122采样点）：$h[114:122]=h[114:122]+\cos (2\pi t\times 0.42)$
  • 第二个突发（134-142采样点）：$h[134:142]=h[134:142]+\cos (2\pi t\times 0.42)$

图2(b) S变换的振幅谱成功检测到了所有四个成分。两个啁啾信号的时频轨迹清晰可见，同时两个高频突发也被准确定位。

图2© STFT（高斯窗）检测到了啁啾信号，但未能检测到两个高频突发。由于时间分辨率不足，两个突发信号的干涉使它们看起来像是一个较低频率的单一信号。

图2(d) Wigner分布在交叉啁啾上显示出非常好的时间和频率分辨率。然而，这种高分辨率的代价是产生了交叉项，这些交叉项严重干扰了图的解释。实际上，交叉项的干涉导致两个突发完全被误表示。

自混叠问题及其解决方案

自混叠的产生机制

在特定$n$处的高斯窗的傅里叶谱称为声音高斯。使用方程(11)，S变换的任何声音都可以计算为时间序列的傅里叶谱$H[m/NT]$与声音高斯的乘积。例如，在方程(11)中，声音高斯为：

$$G(m,n)=e^{-\frac{2{\pi }^2{m}^2}{n^2}}$$

由于时域采样，离散频谱是周期性的。因此，在方程(11)中可以看到，当$n/NT$接近奈奎斯特频率$1/2T$时，声音高斯会重叠到$H[m/NT]$的负频率中。我们称之为自混叠。即使采样率满足奈奎斯特准则，这种自混叠也会发生。声音高斯在高频处变得相当宽，自混叠会引入误差。

特殊奈奎斯特频率

为了最小化自混叠的影响，我们为S变换定义了一个特殊的奈奎斯特频率$s_N$。当声音高斯以$s_N$为中心时，$2\sigma$点位于傅里叶奈奎斯特频率$f_N$。因此：

$$s_n+2\frac{s_n}{2\pi }=f_N$$

最小采样率$T_s$由下式给出：

$$T_s=\frac{\pi +1}{2\pi f_N}$$

与傅里叶奈奎斯特采样率$T_f=1/(2f_N)$相比，这个要求稍微严格一些。当然，在处理实值时间序列时，可以使用时间序列的解析信号进行局部化。这样做的结果是将所有负频率设置为零，从而不存在自混叠。

讨论与应用前景

图1和图2展示了S变换适用的时间序列类别，它们突出了这种方法相对于其他技术的优势。

在图1中，使用包含高频、低频和高频突发的合成时间序列来比较S变换与STFT的性能。S变换的频率依赖分辨率允许检测高频突发，并在长周期信号上显示良好的频率分辨率。STFT的恒定窗口宽度导致所有三个成分的定义较差。STFT中的长周期信号存在自混叠，导致其振幅变化。使用STFT时，在检测短暂的高频信号和低频信号之间存在权衡。

在图2中，将两个高频突发添加到交叉啁啾中。S变换检测到了所有四个成分。STFT可以看到啁啾，但未检测到两个高频成分。STFT没有足够的时间分辨率来解析这两个信号。实际上，两个突发的干涉使其看起来像是较低频率的单个信号。Wigner分布具有非常好的时间和频率分辨率，这在交叉啁啾上很明显。这种分辨率的代价是存在交叉项，阻碍了对此类图的解释。实际上，交叉项的干涉导致两个突发完全被误表示。

结论

从上述例子中可以清楚地看出S变换的优势和局限性。局部化高斯窗的频率倒数依赖性是对STFT中使用的固定宽度窗口的改进。S变换的相位参考到时间原点，提供了关于频谱的有用补充信息，这是CWT中局部参考相位信息所不具备的。然而，与傅里叶变换相关，S变换也受到影响所有离散采样变换的一些相同缺陷的困扰。我们相信S变换的进一步发展将在广泛的学科中找到应用。

附录：S变换的数学推导细节

A. S变换作为傅里叶谱的操作

S变换可以写成对$h(t)$的傅里叶谱$H(f)$的操作。从方程(4)开始：

$$S(\tau ,f)={\int }_{-\infty }^{\infty }h(t)\frac{|f|}{\sqrt{2\pi }}{e}^{-\frac{(t-\tau {)}^2{f}^2}{2}}{e}^{i2\pi ft}dt$$

令$u=t-\tau$，则：

$$S(\tau ,f)={\int }_{-\infty }^{\infty }h(u+\tau )\frac{|f|}{\sqrt{2\pi }}{e}^{-\frac{u^2{f}^2}{2}}{e}^{i2\pi f(u+\tau )}du$$

$$=e^{i2\pi f\tau }{\int }_{-\infty }^{\infty }h(u+\tau )\frac{|f|}{\sqrt{2\pi }}{e}^{-\frac{u^2{f}^2}{2}}{e}^{i2\pi fu}du$$

使用傅里叶变换的卷积定理，我们可以将其写为：

$$S(\tau ,f)={\int }_{-\infty }^{\infty }H(\alpha +f)e^{-\frac{2{\pi }^2{\alpha }^2}{f^2}}{e}^{i2\pi \alpha \tau }d\alpha f\ne 0$$

这个方程(9)显示了S变换如何通过在频率域中使用高斯窗对傅里叶谱进行局部化。

B. 离散S变换的详细推导

从连续S变换的频域表示开始：

$$S(\tau ,f)={\int }_{-\infty }^{\infty }H(\alpha +f)e^{-\frac{2{\pi }^2{\alpha }^2}{f^2}}{e}^{i2\pi \alpha \tau }d\alpha$$

在离散情况下，我们有：

• $\tau \to jT$（时间离散化）
• $f\to n/NT$（频率离散化）
• $\alpha \to m/NT$（积分变量离散化）

代入这些离散化参数：

$$S\left[jT,\frac{n}{NT}\right]=\sum _{m=-\infty }^{\infty }H\left[\frac{m+n}{NT}\right]e^{-\frac{2{\pi }^2{m}^2}{n^2}}{e}^{\frac{i2\pi mj}{N}}$$

由于离散傅里叶变换的周期性，$H[k/NT]$的周期为$N$，因此求和可以限制在一个周期内：

$$S\left[jT,\frac{n}{NT}\right]=\sum _{m=0}^{N-1}H\left[\frac{m+n}{NT}\right]e^{-\frac{2{\pi }^2{m}^2}{n^2}}{e}^{\frac{i2\pi mj}{N}},n\ne 0$$

C. 逆变换的证明

要证明逆变换的正确性，我们需要验证：

$$\sum _{j=0}^{N-1}S\left[jT,\frac{n}{NT}\right]=N\cdot H\left[\frac{n}{NT}\right]$$

从S变换的定义开始：

$$\sum _{j=0}^{N-1}S\left[jT,\frac{n}{NT}\right]=\sum _{j=0}^{N-1}\sum _{m=0}^{N-1}H\left[\frac{m+n}{NT}\right]e^{-\frac{2{\pi }^2{m}^2}{n^2}}{e}^{\frac{i2\pi mj}{N}}$$

交换求和顺序：

$$=\sum _{m=0}^{N-1}H\left[\frac{m+n}{NT}\right]e^{-\frac{2{\pi }^2{m}^2}{n^2}}\sum _{j=0}^{N-1}{e}^{\frac{i2\pi mj}{N}}$$

内部求和是一个几何级数，当$m=0$时等于$N$，否则等于0：

$$\sum _{j=0}^{N-1}{e}^{\frac{i2\pi mj}{N}}=\{\begin{array}{ll}N& \text{if }m=0\\ 0& \text{otherwise}\end{array}$$

因此：

$$\sum _{j=0}^{N-1}S\left[jT,\frac{n}{NT}\right]=N\cdot H\left[\frac{n}{NT}\right]e^0=N\cdot H\left[\frac{n}{NT}\right]$$

这证明了时间平均确实给出了傅里叶谱，从而保证了逆变换的正确性。

D. 声音高斯的频谱分析

在频率$n/NT$处的声音高斯在频域中的表达式为：

$$G_n(m)=e^{-\frac{2{\pi }^2{m}^2}{n^2}}$$

这是一个以零频率为中心、标准差为$\sigma =n/(2\pi )$的高斯函数。随着频率$n$的增加：

• 时域中的高斯窗变窄（宽度$\propto 1/n$）
• 频域中的声音高斯变宽（宽度$\propto n$）

这种倒数关系是S变换频率依赖分辨率的核心。在低频时，S变换使用宽时间窗口（良好的频率分辨率），在高频时使用窄时间窗口（良好的时间分辨率）。

E. 计算复杂度分析

离散S变换的直接计算需要$O(N^2)$次复数乘法运算（对每个$(j,n)$对）。然而，使用快速傅里叶变换（FFT）和卷积定理，可以将计算复杂度降低到$O(N^2\log N)$：

1. 首先计算输入信号的FFT：$O(N\log N)$
2. 对每个频率$n$：
   • 构造声音高斯：$O(N)$
   • 元素乘法：$O(N)$
   • 逆FFT：$O(N\log N)$
3. 总复杂度：$O(N)\times O(N\log N)=O(N^2\log N)$

这使得S变换在实际应用中是可行的，即使对于相对较长的信号也是如此。