抱歉,您的浏览器无法访问本站
本页面需要浏览器支持(启用)JavaScript
了解详情 >

普通傅里叶变换是酉变换,信息完全保存,逆变换可精确重建原信号。由此常引出一个疑问:既然信息一比特未丢,为何仍存在”测不准”?本文的结论是,测不准原理约束的不是信息总量,而是信息的局部可读性。须强调的是,酉性只属于普通傅里叶变换——其加窗版本(短时傅里叶变换)与小波变换都只是等距算子而非酉算子:保范但不满射,而值域的再生核结构恰好给出了测不准原理在时频域最精确的表述。文中先给出 Heisenberg–Weyl 不等式及其证明,继而从相位斜率、幅值谱时序不变性与值域再生核三条机制说明可逆性为何不提供局部可读性,最后以一个跨频段双任务算例验证小波变换的优势边界,并导出一个可直接计算的方法选择判据。

1. 问题的提出:可逆性与可读性

需要区分的是两个不同层次的命题(以下均针对普通傅里叶变换,加窗情形见 4.3 节):

命题 是否成立 所属层次
傅里叶变换可逆,信息不丢失 成立 可逆性
由 $\hat f$ 可精确重建 $f(t)$,任意时刻状态可知 成立 可逆性
由 $f(t)$ 可精确算出 $\hat f(\omega)$,频率成分可知 成立 可逆性
存在联合密度 $P(t,f)$,可读出”$t$ 时刻的频率” 不成立 可读性
可同时把时刻定位到任意精度、把频率分辨到任意精度 不成立 可读性

前三行属可逆性范畴,后两行属可读性范畴,二者互不冲突。常见的推理失误在于把第 2、3 行合取后误推出第 4 行。$f(t)$ 与 $\hat f(\omega)$ 是同一信息的两种互斥表示,不能同时以高分辨率呈现于同一张图上。

2. Heisenberg–Weyl 不等式

2.1 定理表述

设 $f \in L^2(\mathbb{R})$,$\lVert f \rVert_2 = 1$,其傅里叶变换为 $\hat f(\omega)$。定义时间与频率的中心及展宽:

$$ t_0 = \int t \lvert f(t)\rvert^2 dt, \qquad \Delta t^2 = \int (t-t_0)^2 \lvert f(t)\rvert^2 dt $$

$$ \omega_0 = \frac{1}{2\pi}\int \omega \lvert \hat f(\omega)\rvert^2 d\omega, \qquad \Delta\omega^2 = \frac{1}{2\pi}\int (\omega-\omega_0)^2 \lvert \hat f(\omega)\rvert^2 d\omega $$

则恒有

$$ \Delta t \cdot \Delta \omega \geq \frac{1}{2} \qquad\Longleftrightarrow\qquad \Delta t \cdot \Delta f \geq \frac{1}{4\pi} \approx 0.0796 $$

该式在信号分析中称为 Gabor 极限,与量子力学中的位置–动量测不准关系是同一条数学定理,仅替换了共轭变量对。

2.2 证明

不失一般性取 $t_0 = \omega_0 = 0$。由 Parseval 定理及 $\mathcal{F}[f^{\prime}] (\omega) = i\omega\hat f(\omega)$,有 $\Delta\omega^2 = \lVert f^{\prime} \rVert^2$。于是

$$ \Delta t^2 \Delta\omega^2 = \lVert tf \rVert^2 \lVert f^{\prime} \rVert^2 \geq \lvert \langle tf, f^{\prime} \rangle \rvert^2 \geq \lvert \operatorname{Re}\langle tf, f^{\prime} \rangle \rvert^2 $$

第一个不等号为 Cauchy–Schwarz 不等式。注意到 $2\operatorname{Re}\big(\overline{f}f^{\prime}\big) = (\lvert f\rvert^2)^{\prime}$,对实部项分部积分:

$$ \operatorname{Re}\langle tf, f^{\prime} \rangle = \frac{1}{2}\int t \big(\lvert f\rvert^2\big)^{\prime} dt = -\frac{1}{2}\int \lvert f\rvert^2 dt = -\frac{1}{2} $$

代回即得 $\Delta t^2 \Delta\omega^2 \geq 1/4$,证毕。

2.3 等号成立条件

Cauchy–Schwarz 取等要求 $f^{\prime} = -\lambda t f$($\lambda > 0$),解得

$$ f(t) = C e^{-\lambda t^2/2} $$

当且仅当 $f$ 为高斯函数(或其调制平移)时达到下界。取 $f(t) = (\pi\sigma^2)^{-1/4}e^{-t^2/2\sigma^2}$,则 $\hat f(\omega) \propto e^{-\sigma^2\omega^2/2}$,可直接验证:

  • $\lvert f \rvert^2$ 的方差为 $\sigma^2/2$,故 $\Delta t = \sigma/\sqrt2$
  • $\lvert \hat f \rvert^2$ 的方差为 $1/(2\sigma^2)$,故 $\Delta\omega = 1/(\sigma\sqrt2)$
  • 乘积恰为 $1/2$

参数 $\sigma$ 只能在时间与频率之间转移不确定度,不能减少总量。这是 Gabor 变换选用高斯窗的理论依据。

2.4 极端情形与推论

信号 $\Delta t$ $\Delta\omega$
$\delta(t-t_0)$ $0$ $\infty$
$e^{i\omega_0 t}$ $\infty$ $0$
高斯 $e^{-t^2/2\sigma^2}$ $\sigma/\sqrt2$ $1/(\sigma\sqrt2)$

由 Paley–Wiener 定理可进一步推出:不存在时域与频域同时紧支撑的非零函数。这解释了有限长滤波器必然存在过渡带的原因。

3. 时频平面的铺砌

上述约束的几何解释为:时频平面上每个”分辨率单元”的面积存在下界,不同变换只是以不同形状铺砌同一平面。

时频平面的三种铺砌方式

  • 时域采样:单元退化为全频段竖条,$\Delta t \to 0$,$\Delta f \to \infty$
  • STFT / Gabor:单元为形状固定的矩形,全频段共用同一窗长
  • 小波变换:单元形状随尺度伸缩。尺度 $a$ 下小波原子为 $a^{-1/2}\psi\left(\frac{t-b}{a}\right)$,时间展宽正比于 $a$,频率展宽反比于 $a$,乘积与 $a$ 无关

三者的单元面积同受 $\Delta t\cdot\Delta f \geq 1/(4\pi)$ 约束。小波并未突破该约束,只是把分辨率按频率重新分配。该下界在加窗分析中的精确来源是变换值域的再生核,即分析窗自身的模糊函数,详见 4.3 节。

4. 可逆性为何不提供局部可读性

4.1 机制一:时间信息编码于相位的斜率

考察最极端的时域局部事件——瞬时冲击 $\delta(t-t_0)$,其傅里叶变换为

$$ \mathcal{F}[\delta(t-t_0)] (\omega) = e^{-i\omega t_0} $$

幅值谱恒为 1,在整条频率轴上完全平坦。”事件发生于 $t_0$”这条信息一点未丢,但它不驻留于任何单个频率分量,而是以相位关于频率的斜率形式存在:

$$ \frac{d\varphi}{d\omega} = -t_0 $$

测量一条直线的斜率需要在自变量上有足够长的基线。若观测带宽仅为 $B$,则

$$ \delta t_0 \sim \frac{1}{B} $$

这正是测不准原理的另一副面孔:时域定位精度由频域观测跨度决定

4.2 机制二:幅值谱无法区分时序

构造一对对照信号:$A(t)$ 为先低频后高频的调制信号,$B(t) = A(t_{\rm end} - t)$ 为其时间反转。由傅里叶变换的时间反转性质,对实信号有

$$ \hat B(\omega) = \overline{\hat A(\omega)} e^{-i\omega t_{\rm end}} \qquad\Longrightarrow\qquad \lvert \hat B(\omega) \rvert \equiv \lvert \hat A(\omega) \rvert $$

两者的幅值谱逐点精确相等。无论频率轴细分到何种程度,幅值谱都无法判定事件的先后次序,差异 100% 存放于相位之中。逆变换固然能把二者完美区分,但那要求使用全频段、全精度的完整相位做一次全局重建;而”看一眼频谱即读出时刻”属于局部读取,恰是被禁止的操作。

4.3 机制三:加窗分析的值域受再生核约束

前两小节讨论的是普通傅里叶变换。转入加窗分析时算子性质发生改变,而这一改变给出了测不准原理在时频域最精确的表述。

设窗 $g$ 满足 $\lVert g\rVert_2 = 1$,短时傅里叶变换(取 $e^{-2\pi i t\omega}$ 约定)为

$$ V_g f(b, \omega) = \int f(t) \overline{g(t-b)} e^{-2\pi i t\omega} dt $$

Moyal 正交关系

$$ \langle V_{g_1} f_1, V_{g_2} f_2\rangle_{L^2(\mathbb{R}^2)} = \langle f_1, f_2\rangle \overline{\langle g_1, g_2\rangle} $$

取 $g_1 = g_2 = g$、$f_1 = f_2 = f$ 即得 $\lVert V_g f\rVert_{L^2(\mathbb{R}^2)} = \lVert f\rVert_{L^2(\mathbb{R})}$,故 $V_g$ 保范。但 $V_g$ 把一元函数映为二元函数,值域必为 $L^2(\mathbb{R}^2)$ 的真闭子空间,不可能满射。于是

$$ V_g^{\ast}V_g = I, \qquad V_g V_g^{\ast} = P \neq I $$

其中 $P$ 为到值域的正交投影。$V_g$ 是等距算子而非酉算子——酉要求双射等距。连续小波变换 $W_\psi$ 情形相同。

变换 保范 满射
傅里叶变换 $\mathcal{F}$
短时傅里叶变换 $V_g$
连续小波变换 $W_\psi$

值域的再生核。 值域是一个再生核 Hilbert 空间:任何合法的时频图 $F = V_g f$ 必满足不动点方程

$$ F(z) = \int_{\mathbb{R}^2} F(z^{\prime}) \langle \pi(z^{\prime})g, \pi(z)g\rangle dz^{\prime} $$

其中 $\pi(z)$ 为时频平移算子,核的模 $\lvert\langle \pi(z^{\prime})g, \pi(z)g\rangle\rvert = \lvert V_g g(z-z^{\prime})\rvert$ 恰为窗自身的模糊函数(ambiguity function)。由此得到测不准原理在 STFT 域的精确表述:

任何时频图都是自身与窗的模糊函数的卷积不动点,因而其集中度不可能细于窗自身的时频足迹。

对高斯窗,该足迹的面积恰为下界 $1/(4\pi)$。画不出任意锐利的时频图,并非算法不足,而是锐利的图根本不在值域之内。

数值核验。 取 $N = 64$、高斯窗、跳步为 1 的全密度离散 Gabor 系统($V$ 把 $N$ 维向量映为 $N\times N$ 矩阵),结果如下(代码见附录 A.1):

检验项 结果 含义
$\lVert Vf\rVert^2/\lVert f\rVert^2$ 恒为 $64 = N\lVert g\rVert^2$,与 $f$ 无关 保范(等距)
$\lVert V^{\ast}Vf/N - f\rVert/\lVert f\rVert$ $4.8\times10^{-15}$ 左逆存在,可精确重建
随机取 $R$,$\lVert PR-R\rVert/\lVert R\rVert$ $0.991$ $R$ 几乎完全不在值域内
$\lVert P(PR)-PR\rVert/\lVert PR\rVert$ $5.6\times10^{-15}$ $P$ 幂等,确为投影
合法时频图 $F$:$\lVert PF-F\rVert/\lVert F\rVert$ $4.8\times10^{-15}$ 再生核性质成立
值域维数 / 环境空间维数 $64/4096 = 1.56%$ 值域是极薄的子空间

即:时频平面上 98% 以上的函数不是任何信号的短时傅里叶变换

离散化后的 STFT 构成 Gabor 框架,一般是冗余的;在临界密度处,Balian–Low 定理(文献 7)进一步表明良好的时频局部化与 Riesz 基不可兼得。

4.4 一个直观类比

设某路口信号灯周期恒为 60 s。连续记录 1 h 后做 FFT,谱线出现于 $1/60 = 16.67\ \mathrm{mHz}$,逆变换可精确恢复任意时刻的红绿状态。此处”时间与频率都可知”成立——原因并非傅里叶变换绕过了测不准原理,而是该信号平稳,频率不随时间变化,联合定位问题从未被提出。

改为自适应信号灯,于某时刻由 60 s 周期切换为 40 s 周期,并提出联合问题:”周期是多少,以及它在何时改变?”区分两周期要求频率分辨率优于

$$ \Delta f = \frac{1}{40} - \frac{1}{60} = 8.33\ \mathrm{mHz} $$

即使采用最优高斯窗,取理论下界 $\Delta t \geq 1/(4\pi\Delta f)$ 得 $\Delta t \geq 9.55\ \mathrm{s}$。“周期由 60 s 变为 40 s,切换发生在某时刻 ±1 s 之内”这一论断在数学上不可能同时成立,与采样率、数据量、算法先进程度均无关。

5. 数值算例:跨频段双任务

5.1 两类约束的区分

讨论小波是否占优之前,须先区分两条性质完全不同的约束:

约束 内容 可否解除
物理约束 单个分析原子必须满足 $\Delta t\cdot\Delta f \geq 1/(4\pi)$;其在时频域的精确形式即 4.3 节的再生核约束 不可解除
人为约束 STFT 在全频段共用同一窗,即所有原子形状相同 可解除

由此得到判定小波能否占优的关键:

  • 两项任务同处一个频段时,二者竞争同一分辨率单元,须由同一个原子承担,触及物理约束。此时相容条件为 $\Delta t_{\rm sep}\cdot\Delta f_{\rm sep} \gtrsim 1/\pi$,不满足则任何方法(含小波)皆无解。
  • 两项任务分处相距甚远的两个频段时,二者由不同尺度的原子分别承担,各原子只需各自满足物理约束,彼此不竞争。此时 STFT 的失败纯粹源于人为约束,小波通过对不同频率使用不同原子将其解除。

本节构造的算例属于后者。

5.2 信号设计

采用旋转机械故障诊断场景($f_s = 10\ \mathrm{kHz}$,时长 2 s):

  • 低频任务——两根转轴的基频 24.0 Hz 与 30.0 Hz($\Delta f_{\rm sep} = 6\ \mathrm{Hz}$)。提问:是哪根轴失衡
  • 高频任务——轴承缺陷激发的 3000 Hz 共振冲击,发生于 $t = 1.000\ \mathrm{s}$ 与 $t = 1.015\ \mathrm{s}$($\Delta t_{\rm sep} = 15\ \mathrm{ms}$)。提问:冲击发生在何时

两项提问格式相同(”是哪一个”与”在何时”),但对分析参数的要求方向相反,且位于相距 111 倍的两个频段

跨频段双任务信号

乘积 $\Delta t_{\rm sep}\cdot\Delta f_{\rm sep} = 0.015 \times 6 = 0.090$,小于单原子双任务阈值 $1/\pi = 0.318$,故任何单一窗长的 STFT 均无解;但大于绝对下界 $1/(4\pi) = 0.0796$,且两项任务分处不同频率,因此并未被物理约束禁止。

5.3 可行域解析

STFT 必须以单一窗长同时满足两项要求。以 Hann 窗主瓣宽度 $2/T_{\rm win}$ 为分辨判据:

$$ T_{\rm win} \geq \frac{2}{\Delta f_{\rm sep}} = 333\ \mathrm{ms} \quad\wedge\quad T_{\rm win} \leq \Delta t_{\rm sep} = 15\ \mathrm{ms} \quad\Longrightarrow\quad \text{可行域} = \varnothing $$

下限为上限的 22 倍。

CWT 的自由参数是品质因数 $Q$(Morlet 小波中 $Q = \omega_0$)。按幅值约定,其分辨率随频率缩放:

$$ \Delta f = \frac{f}{Q}, \qquad \sigma_t = \frac{Q}{2\pi f}, \qquad \sigma_t \cdot \Delta f = \frac{1}{2\pi}\ \text{(与频率无关)} $$

两项要求分别施加于各自所在频率,故互不冲突:

$$ Q \geq \frac{2 f_L}{\Delta f_{\rm sep}} = \frac{2\times 27}{6} = 9.0 \quad\wedge\quad Q \leq \pi f_H \Delta t_{\rm sep} = \pi\times 3000\times 0.015 = 141 $$

$$ \Longrightarrow\quad \text{可行域} = [9,\ 141] \neq \varnothing $$

取 $Q = 24$,实际分辨率为:27 Hz 处 $\Delta f = 1.12\ \mathrm{Hz} \ll 6\ \mathrm{Hz}$;3000 Hz 处 $\sigma_t = 1.27\ \mathrm{ms} \ll 15\ \mathrm{ms}$,两项均有充分裕度。

参数可行域对比

5.4 数值结果

对每种方法给出一张统一的全频段时频图(对数频率轴,时间窗 180 ms),使两项任务在同一张图内同时可判读:青色虚线标出 24 与 30 Hz,绿色虚线标出 1.000 与 1.015 s。

四种方法的统一时频图

方法 分辨 24/30 Hz 分辨 1.000/1.015 s 双任务
STFT 12.8 ms 失败
STFT 102 ms 失败
STFT 819 ms 失败
CWT $Q=24$ 通过

判据采用 Rayleigh 谷深准则:两主峰之间的谷值须低于峰值的 0.7 倍,方判为分辨。该准则排除幅值涟漪造成的伪双峰。

只有 CWT 一列同时具备两项特征:低频端两条清晰谱线分居 24 与 30 Hz,高频端两条清晰竖线分居 1.000 与 1.015 s。三种 STFT 窗长各失败于至少一项,且失败方向随窗长单调迁移——短窗保时间、长窗保频率,中间窗两头落空。

定量切片

定量切片中,CWT 曲线是唯一在左右两图中均呈现双峰的曲线。

5.5 图像细节与方法学要点

其一,819 ms 面板中 3000 Hz 处的贯穿全宽弥散带。 冲击能量被展宽至 819 ms,远超 180 ms 的显示窗,故呈均匀水平带而非竖线。信息并未丢失,只是时间定位被彻底抹平。

其二,102 ms 与 819 ms 面板中高频区的水平细纹。 条纹间距为 $1/\Delta t_{\rm sep} = 66.7\ \mathrm{Hz}$,是两次冲击的干涉条纹。长于冲击间隔的窗把”两次冲击”记录为”一次事件加频域条纹”——同一事实的另一种编码形式。

其三,两种描述均正确。 12.8 ms 窗给出的描述是”单个 27 Hz 分量,幅值以 6 Hz 拍起伏”(拍周期 167 ms,在时域波形图中清晰可见);819 ms 窗给出的描述是”两个独立分量,分别位于 24 与 30 Hz”。二者描述同一信号,均无错误。测不准原理并未销毁信息,只决定信息以何种形式呈现——这与第 4 节的论述是同一现象。

其四,补零不能提高分辨率。 判定计算中三种窗长均补零至 16384 点,频率栅格统一为 $f_s/N_{\rm FFT} = 0.61\ \mathrm{Hz}$,远细于 6 Hz 间隔。12.8 ms 窗的频率切片光滑、采样充分,却依然只有单峰。分辨率由窗长决定,与 FFT 点数无关——补零仅对频谱做 sinc 插值,不引入新信息。

6. 一般性判据:可行域宽度比

设两项任务分别位于频率 $f_L$(要求频率分辨率 $\Delta f_{\rm sep}$)与 $f_H$(要求时间分辨率 $\Delta t_{\rm sep}$),$f_H > f_L$。由 5.3 节两式定义可行域宽度比

$$ R = \frac{Q_{\max}}{Q_{\min}} = \frac{\pi}{2}\cdot\frac{f_H}{f_L}\cdot \Delta t_{\rm sep}\cdot\Delta f_{\rm sep} \approx 1.571\ \frac{f_H}{f_L}\ \Delta t_{\rm sep}\ \Delta f_{\rm sep} $$

判定规则为:

  • $R < 1$:小波无解
  • $R \approx 1$:可行域退化,解不稳健
  • $R \gg 1$:可行域宽裕,小波稳健占优

等价地,小波有解要求

$$ \frac{f_H}{f_L} \geq \frac{2}{\pi}\cdot\frac{1}{\Delta t_{\rm sep}\cdot\Delta f_{\rm sep}} \approx \frac{0.637}{\Delta t_{\rm sep}\cdot\Delta f_{\rm sep}} $$

代入本文算例:$\Delta t_{\rm sep}\cdot\Delta f_{\rm sep} = 0.090$,要求 $f_H/f_L \geq 7.1$,而实际为 111,故 $R = 15.7$,裕度充分。作为对照,若两项任务同处 100 Hz 附近($f_H/f_L \approx 2.9$)且 $\Delta t_{\rm sep}\cdot\Delta f_{\rm sep} = 0.16$,则 $R = 0.74 < 1$,小波同样无解——这正是 5.1 节所述”竞争同一分辨率单元”的情形。

该式表明:两项任务的频率比越大、各自的分辨要求越宽松,小波的可行域越宽。$R$ 可在设计分析方案之前直接由任务指标算出,无需试算。

7. 结论与适用范围

小波的优势不是”突破了测不准原理”,而是”解除了 STFT 全频段共用一个窗这一人为约束”。 据此可给出明确的方法选择准则:

  • 应当用小波:信号具有多尺度结构,低频分量需要高频率分辨率、高频分量需要高时间分辨率,且 $R \gg 1$。典型场景为旋转机械故障诊断、地震初至与面波联合分析、心电信号的 ST 段与 QRS 波分析、瞬变电磁响应。
  • 不必用小波:所有关注对象处于同一频段(此时恒 $Q$ 特性无用武之地,长窗 STFT 更优),或信号本身平稳(直接用 FFT 即可)。
  • 小波亦无解:$R < 1$,即两项任务竞争同一分辨率单元。此时需转向同步压缩、重排谱等非线性方法,或更换传感方式——例如以事件时间戳直接记录,事件记录不做时频分解,因而不受该约束。

附录:可复现代码

环境为 Python 3.14 + NumPy 2.4 + SciPy 1.18。

A.1 算子性质核验(4.3 节)

# -*- coding: utf-8 -*-
"""离散 STFT 是等距算子但非酉算子"""
import numpy as np

rng = np.random.default_rng(7)
N = 64
k = np.arange(N)
E = np.exp(-2j * np.pi * np.outer(k, k) / N)          # E[m, t] = exp(-2i pi m t / N)
g = np.exp(-0.5 * ((k - N / 2) / (N / 8.0)) ** 2)
g = g / np.linalg.norm(g)                              # 单位范数高斯窗
G = np.stack([np.roll(g, j) for j in range(N)], axis=1)   # G[t, n] = g[t-n]


def V(f):
    """离散 STFT:N 维向量 -> N x N 复矩阵(行为频率,列为时移)"""
    return E @ (f[:, None] * np.conj(G))


def Vstar(F):
    """伴随算子"""
    return np.einsum("mn,mt,tn->t", F, np.conj(E), G)


f1 = rng.standard_normal(N) + 1j * rng.standard_normal(N)
f2 = rng.standard_normal(N) + 1j * rng.standard_normal(N)
c = N * np.linalg.norm(g) ** 2

r1 = np.linalg.norm(V(f1)) ** 2 / np.linalg.norm(f1) ** 2
r2 = np.linalg.norm(V(f2)) ** 2 / np.linalg.norm(f2) ** 2
print("保范  : {:.10f}  {:.10f}   理论 N||g||^2 = {:.10f}".format(r1, r2, c))
print("左逆  : ||V*Vf/c - f||/||f|| = {:.3e}".format(
    np.linalg.norm(Vstar(V(f1)) / c - f1) / np.linalg.norm(f1)))

R = rng.standard_normal((N, N)) + 1j * rng.standard_normal((N, N))   # 时频面上任取一函数
PR = V(Vstar(R)) / c
print("非满射: ||PR - R||/||R||     = {:.4f}".format(
    np.linalg.norm(PR - R) / np.linalg.norm(R)))
print("幂等  : ||P(PR) - PR||/||PR|| = {:.3e}".format(
    np.linalg.norm(V(Vstar(PR)) / c - PR) / np.linalg.norm(PR)))

F1 = V(f1)
print("再生核: ||PF - F||/||F||     = {:.3e}".format(
    np.linalg.norm(V(Vstar(F1)) / c - F1) / np.linalg.norm(F1)))
print("值域维数 {} / 环境空间维数 {} = {:.2f}%".format(N, N * N, 100.0 * N / N ** 2))

A.2 跨频段双任务判定(5.4 节)

以下脚本复现正文的可行域解析与判定表,不含绘图部分。

# -*- coding: utf-8 -*-
"""跨频段双任务的数值判定

低频任务:分辨 24 / 30 Hz 两个转轴分量         -> 要求长窗
高频任务:分辨 3000 Hz 处相隔 15 ms 的两次冲击 -> 要求短窗
结论:任何单一 STFT 窗长均失败;Morlet 小波在 Q in [9, 141] 内成功。
"""
import numpy as np
from scipy.signal import ShortTimeFFT
from scipy.signal.windows import hann


def stft_mag(x, fs, N, hop, mfft):
    """Hann 窗 STFT 幅值谱。补零至 mfft 点,使频率栅格远细于分辨率。"""
    S = ShortTimeFFT(hann(N, sym=False), hop=hop, fs=fs, mfft=mfft, scale_to="magnitude")
    return np.abs(S.stft(x)), S.f, S.t(x.size)


def cwt_morlet(x, fs, freqs, Q):
    """Morlet 连续小波变换。幅值保持归一化:单位幅值正弦在其自身尺度上给出 |W| = 1。"""
    Xf = np.fft.fft(x)
    w = 2 * np.pi * np.fft.fftfreq(x.size, d=1 / fs)
    pos = w > 0
    out = np.empty((freqs.size, x.size))
    for i, fc in enumerate(freqs):
        a = Q / (2 * np.pi * fc)
        out[i] = np.abs(np.fft.ifft(Xf * (2.0 * np.exp(-(a * w - Q) ** 2 / 2) * pos)))
    return out


def resolved(y, thr=0.5, dip=0.7):
    """Rayleigh 谷深判据:两主峰之间谷值须低于峰值的 dip 倍,方判为分辨。"""
    idx = [i for i in range(1, len(y) - 1)
           if y[i] > y[i - 1] and y[i] >= y[i + 1] and y[i] > thr]
    if len(idx) < 2:
        return False
    a, b = sorted(sorted(idx, key=lambda i: -y[i])[:2])
    return y[a:b + 1].min() <= dip * min(y[a], y[b])


def nrm(v):
    return v / v.max()


# ---------------------------------------------------------------- 信号构造
fs = 10000.0
t = np.arange(0, 2.0, 1 / fs)
FL1, FL2, FH = 24.0, 30.0, 3000.0      # 两转轴基频;轴承共振频率
TH1, TH2 = 1.000, 1.015                # 两次冲击时刻
DF, DT = FL2 - FL1, TH2 - TH1          # 6 Hz;15 ms
FL = 0.5 * (FL1 + FL2)

x = np.sin(2 * np.pi * FL1 * t) + np.sin(2 * np.pi * FL2 * t)
for t0 in (TH1, TH2):
    m = t >= t0
    x[m] += 2.5 * np.exp(-(t[m] - t0) / 6e-4) * np.sin(2 * np.pi * FH * (t[m] - t0))
x += 0.01 * np.random.default_rng(1).standard_normal(t.size)

# ---------------------------------------------------------------- 可行域解析
T_LO = 2.0 / DF                 # STFT 低频任务:Hann 窗主瓣 2/T_win <= df_sep
T_HI = DT                       # STFT 高频任务:T_win <= dt_sep
Q_MIN = 2 * FL / DF             # CWT 低频任务:2*Df = 2*f_L/Q <= df_sep
Q_MAX = np.pi * FH * DT         # CWT 高频任务:2*sigma_t = 2*Q/(2*pi*f_H) <= dt_sep
Q_USE = 24.0

print("频率比 f_H/f_L = {:.0f},  dt_sep * df_sep = {:.3f}".format(FH / FL, DT * DF))
print("STFT: T_win >= {:.0f} ms 且 T_win <= {:.0f} ms  ->  空集(相差 {:.0f} 倍)"
      .format(T_LO * 1e3, T_HI * 1e3, T_LO / T_HI))
print("CWT : Q >= {:.1f} 且 Q <= {:.1f}  ->  [{:.0f}, {:.0f}],R = {:.1f},取 Q = {:.0f}"
      .format(Q_MIN, Q_MAX, Q_MIN, Q_MAX, Q_MAX / Q_MIN, Q_USE))
print("Q = {:.0f} 时:{:.0f} Hz 处 Df = {:.2f} Hz;  {:.0f} Hz 处 sigma_t = {:.2f} ms\n"
      .format(Q_USE, FL, FL / Q_USE, FH, Q_USE / (2 * np.pi * FH) * 1e3))

# ---------------------------------------------------------------- 数值判定
LO, HI = (15.0, 40.0), (2000.0, 4000.0)
rows = []
for N in (128, 1024, 8192):
    S, f, tt = stft_mag(x, fs, N, hop=20, mfft=16384)
    fm = (f >= LO[0]) & (f <= LO[1])
    r1 = resolved(nrm(S[fm, int(np.argmin(np.abs(tt - 0.60)))]))
    bm = (f >= HI[0]) & (f <= HI[1])
    env = nrm(np.sqrt((S[bm] ** 2).mean(axis=0)))
    tm = (tt >= 0.97) & (tt <= 1.05)
    rows.append(("STFT {:.0f} ms".format(N / fs * 1e3), r1, resolved(env[tm])))
    del S

freqs = np.geomspace(15.0, 4200.0, 340)
W = cwt_morlet(x, fs, freqs, Q=Q_USE)
lom = (freqs >= LO[0]) & (freqs <= LO[1])
him = (freqs >= HI[0]) & (freqs <= HI[1])
r1 = resolved(nrm(W[lom, int(np.argmin(np.abs(t - 0.60)))]))
env = nrm(np.sqrt((W[him] ** 2).mean(axis=0)))
tm = (t >= 0.97) & (t <= 1.05)
rows.append(("CWT Q={:.0f}".format(Q_USE), r1, resolved(env[tm])))

print("{:<14s}{:>10s}{:>10s}{:>9s}".format("方法", "低频任务", "高频任务", "双任务"))
print("-" * 46)
for name, r1, r2 in rows:
    print("{:<14s}{:>10s}{:>10s}{:>9s}".format(
        name, "分辨" if r1 else "合并", "分辨" if r2 else "合并",
        "通过" if (r1 and r2) else "失败"))
print("-" * 46)
print("通过双任务的方法数:{}".format(sum(1 for _, a, b in rows if a and b)))

参考文献

  1. Gabor D. Theory of communication. Journal of the Institution of Electrical Engineers, 1946, 93(26): 429–457.
  2. Folland G B, Sitaram A. The uncertainty principle: a mathematical survey. Journal of Fourier Analysis and Applications, 1997, 3(3): 207–238.
  3. Daubechies I. Ten Lectures on Wavelets. Philadelphia: SIAM, 1992.
  4. Mallat S. A Wavelet Tour of Signal Processing: The Sparse Way. 3rd ed. Burlington: Academic Press, 2008.
  5. Torrence C, Compo G P. A practical guide to wavelet analysis. Bulletin of the American Meteorological Society, 1998, 79(1): 61–78.
  6. Donoho D L, Stark P B. Uncertainty principles and signal recovery. SIAM Journal on Applied Mathematics, 1989, 49(3): 906–931.
  7. Benedetto J J, Heil C, Walnut D F. Differentiation and the Balian–Low theorem. Journal of Fourier Analysis and Applications, 1995, 1(4): 355–402.
  8. Daubechies I, Lu J, Wu H-T. Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool. Applied and Computational Harmonic Analysis, 2011, 30(2): 243–261.
  9. 抖音视频. https://v.douyin.com/YxlCjq6-NnI/

评论