当前位置：首页 > java >正文

信号(瞬时)频率求解与仿真实践(1)

java 2025/6/25 23:15:52

引言

信号频率求解的重要性无需多言，其在各方各面都有迫切的需求。比如在雷达信号处理领域：FMCW体制雷达去斜接收下，我们需要对各脉冲做FFT处理得到频率信息才能获取目标的距离信息；我们需要在慢时间维做FFT处理得到多普勒频率才能获取目标的速度信息。

最近在课题的学习中碰到了有关信号分解的话题，做了些调研和仿真，想着把该话题系统地梳理一下，于是形成了本博文。我预设了两篇文章探讨该话题，本文是第一篇，主要是搭建框架以及对瞬时频率分析方法做介绍，第二篇博文(对应文中3.2节)将主要讨论联合时频分析法。

这两篇博文的内容并没有涉及具体的课题应用，理论和分析都很简单。如有错误之处欢迎读者留言指正。

Blog

2025.6.6 本文第一次撰写

一、信号分类及各类别信号处理方法概述

本文针对的是信号频率的求解，从频率的角度出发，我们可以将信号分成：频率时变信号和频率时不变信号。顾名思义，频率时变信号就是指信号的频率在信号的采集时间段内是在随时间变化的，最简单的如chirp信号,其频率f(t)为：

（1-1）

f0 是信号采集时的起始频率，Kr是调频斜率。其对应的时域表达式为（如果是双路采样）：

（1-2）

频率时不变信号是指信号频率至少在整个采样时间段内是不随时间变化的，最简单的如单频正弦信号：

(1-3)

该信号在整个采样时间区间内频率值都是f0。

对于频率时不变信号，瞬时频率这个概念没什么意义，因为频率每时每刻都一样。此类信号求其频率时，我们可以直接通过傅里叶变换获取信号频谱，并进而从频谱中得到信号的频率值(且不管信号中有多少个频率分量，只要对信号的采样率和采样时间设置合理(涉及到最大可见频率值以及频率分辨率)，则都能看到这些频率分量形成的波峰)。

对于频率时变信号，则不能使用传统的傅里叶变换求解，因为传统的傅里叶变换是没有时间维度信息的。这种情况下一般有两种处理方法：1、瞬时频率法； 2、联合时频分析法。

瞬时频率法，就是求解信号相位对时间的导数来获取瞬时频率。对于单分量信号(信号中每个时间点下只有一个频率分量)，如果采集到的是复数值，那么可以直接求相位然后求对时间的导数即可；如果采集到的是实数值，则可以使用希尔伯特变换构造解析信号(可以简单理解成实数信号对应的复数信号)，然后求该解析信号的相位，进而对时间求导得到频率信息。但是这种方法无法求解当信号中有多个频率分量的情况（因为当信号中有多个频率分量时，它们会叠加在一起，直接基于相位求导获取的频率信息是多个分量混叠的结果），如果信号中有多个频率分量，此时需要先将信号分解成多个单分量的信号，再分别对各个单分量信号进行希尔伯特变换，求解其瞬时频率，最终叠加各个分量的处理结果，获得该多分量信号的时频分布。典型的信号分解方法如：经验模态分解(Empirical Mode Decomposition, EMD)。【事实上，先做信号分解，后使用希尔伯特法得到各分量的解析信号，并求各解析信号的瞬时频率，该流程论述的就是所谓的希尔伯特-黄变换(Hilbert-Huang Transform，HHT)】

联合时频分析法是另一种解决思路，它不需要求解信号的相位以及做多分量信号的分解。而是同时在时间和频率维度处理信号，经典的方法如：短时傅里叶变换(STFT，线性时频变换的一种)、小波变换、CVD谱、Cepstrogram谱、以及使用不同核函数解决信号自相关交叉项干扰的各类非线性时频变换法，如：维格纳 - 维利分布(Wigner-Ville distribution，WVD)等等。 联合时频分析法可以直接给出信号在二维的时间-频率域的能量分布，不仅可用于单分量信号也可以直接用于多分量信号的时频分析。这部分内容更具体的我将在下一篇文章中展开。

把前面的内容用一张图概括，大概如下图所示：

图1.1 信号(瞬时)频率分析与处理框架

一些更细节的理论内容会在后文与仿真的结果一起展开。

二、频率时不变信号频率求解

求解方法：直接FFT处理得到频谱，（不管信号中有多少个频率分量）。至于采样率是否满足最大频率要求、采样时间是否满足最小频率分辨率的要求等，不再本文的讨论范围。

本章的仿真中，我预设了两个频率分量：20Hz和40Hz，采样率为60Hz，双路采样，采样点数为512个点，加SNR为20dB的白噪声，得到的时域和频域结果如下：

图2.1 信号时域

图2.2 信号频域

三、频率时变信号频率求解

3.1 瞬时频率分析法

瞬时频率法就是通过求解信号的相位，然后用相位对时间求导来得到信号的瞬时频率。最核心的在于获取信号的相位，而要求解信号的相位只能是复数信号。此外，还需要注意该方法只能处理单分量的信号，对于多分量信号我们需要先进行信号分解。

3.1.1 单分量信号

单分量信号下不用考虑多种频率信号的混叠，但是要分成复数信号和实数信号两种情况讨论。

A.采集的是复数信号

这是最简单的情况，直接求解信号的相位并求其对时间的导数即可。

本节的仿真中，我构造了一个chirp信号，其起始频率为20Hz,调频斜率为50Hz/s，设置采样率400Hz，采集512个点(双路采样)，并加SNR为20dB的白噪声，仿真的结果如下：

图3.1 信号时域

如果直接进行FFT处理，我们能看到信号所包含的频率信息，但是无法获知频率在时间维度的分布：

图3.2 直接FFT处理

(无法获知频率的时间分布)

求该信号的相位(需要对相位解跳变，信号相位是在-pi至pi之间变化的，为了看到相位随时间真实的变化信息，需要进行unwarp解跳变操作)：

图3.3 信号相位求解结果(解跳变前后)

随后将解跳变后的相位对时间求导数，并除以(2*pi)得到信号频率随时间的变化：

图3.4 信号的瞬时频率

毛刺是由噪声对相位的扰动产生的。结果符合预设条件。

B.采集的是实数信号

我们无法直接从实数信号中获取相位信息。一种典型的将实数信号转成其相对应的复数信号的方法是希尔伯特变换。

希尔伯特变换就是一个数学工具，可以帮助我们把实数信号转换为解析信号(复数形式)，从而可以提取其瞬时幅度、相位以及频率信息。该变换定义如下：对于实数信号x(t),对其做希尔伯特变换，得到:

（3-1）

式中，p.v.为柯西主值（避免积分奇点）。从时域上来看，对信号做希尔伯特变换就是将之与函数做卷积。得到的其实就是x(t)的正交分量(相位延迟90°后的信号)。随后我们可以得到解析信号为：

（3-2）

这就是复数信号了，此时我们可以用前面A中的方法，求相位、求相位对时间的导数就可以得到瞬时频率信息。

至于如何实现，Matlab有函数：hilbert，输入实数信号输出的就是其对应的解析信号。 或者也可以自编代码实现： 将原始信号做FFT，然后将其正频率部分乘以2，其余置0，再IFFT回时域也可以近似得到其解析信号(考虑到1/(pi*t)在t=0处存在奇异点，以及我们无法在无限长的时间内进行卷积，所以不建议直接在时域上做卷积处理。上面给出的是频域处理方法：将实数信号的FFT结果中负频率部分置0后频域失去共轭对称性，再IFFT回去则可以得到复数信号，之所以将正频率部分乘以2是为了保证能量守恒？解析信号中引入了希尔伯特变换后的值，我们需要给正频率分量的幅度加倍才能保证转换前后的信号能量守恒(加上了希尔伯特变换值)。有关希尔伯特变换更细节的理论说明可以参考其它资料，这不是本文的重点。

本节的仿真参数与前面A中采集的是复数信号参数一致，不过这里采集单路信号。仿真的结果如下：

图3.5 希尔伯特变换前后的信号时域

从上图最下面的就是信号希尔伯特变换的结果，可以看到其本质就是原始信号相位延迟90°后的结果。

随后，与前面复数信号的处理流程一样，获取解析信号的相位、并求解相位对时间的导数得到信号的瞬时频率信息：

图3.6 解析信号相位解跳变前后

图3.7 信号频率随时间的变化

结果基本符合预设条件。

3.1.2 多分量信号

应对多分量信号，最直接的想法就是把各个分量单独拎出来(信号分解)。其实傅里叶变换的本质就是信号分解，但它最致命的缺点是无法得到频率的时间信息(我们可以从FFT的结果中看到全部的频率信息，但是不知道这些频率在时间维度是怎么样分布的)。短时傅里叶变换可以解决我们对于观察信号频率在时间上分布的需求，但它是以牺牲频率分辨率为代价的。再后来有小波变换，可是小波变换需要选用特定的小波基，一旦确定了小波基，在整个分析过程中将无法更换，即使该小波基在全局可能是最佳的，但在某些局部可能并不是，所以小波分析的基函数缺乏适应性[1]。

如果多分量信号得到有效分解，那么后续就针对各个分量实践前面3.1.1的处理流程即可，所以本节的重点在于介绍如何分解，特别地，介绍EMD。

EMD是什么？怎么实现？

经验模态分解（Empirical Mode Decomposition, EMD）是一种自适应信号处理方法，由黄锷（Norden E. Huang）于1998年提出，适用于分析非线性、非平稳信号。其核心思想是将复杂信号分解为有限个本征模态函数（Intrinsic Mode Functions, IMF）（不要被名字误导，其实就是信号分量），每个IMF代表信号的一个固有振荡模式。

不过本征模态函数分量有两个约束条件：

1. 在整个数据段内，你分解出来的本征模态函数的极值点个数和过零点的个数必须相等或相差最多不能超过一个。

2. 任意时刻，由局部极大值点形成的上包络线和由局部极小值点形成的下包络线的平均值为零，即上、下包络线相对于时间轴局部对称。

怎么把该复杂的多分量信号分解成多个本征模态函数分量？(下面的内容主要来自参考资料[1])：

1）根据原始信号上下极值点，分别画出上、下包络线。

2) 求上、下包络线的均值，得到均值包络线。

3）原始信号减均值包络线，得到中间信号(residue)。

4）判断该中间信号是否满足IMF的两个条件，如果满足，该信号就是一个IMF分量；如果不是，以该信号为基础，重新做1）~ 4）的分析。IMF分量的获取通常需要若干次的迭代。

5）使用上述方法得到第一个IMF后，用原始信号减IMF1，作为新的原始信号，再通过1）~4）的分析，可以得到IMF2，以此类推，完成EMD分解。（上面的图片内容来自参考资料[2]，读者可以通过链接下载，该PPT描述了一个完整的信号EMD分解过程，我也一并附在仿真代码链接中了）

EMD最终会输出多少个IMF？

这由EMD算法的迭代截止条件决定，比如：1、你可以设置：如果迭代过程中其上下包络的局部均值接近零（如均值绝对值<某个阈值）则终止迭代，留下的就用残差项来代替； 2、你可以自行预设最大的IMF个数这个超参，当计算得到预设的最大IMF个数后不再迭代； 3、或者残余项变成常量；需要注意的是，过少的IMF可能导致模态混叠（不同频率成分未分离）；过多的IMF可能引入虚假分量（如噪声被分解为多个IMF）。至于到底要多少个，结合实际应用。

EMD有哪些应用

1、信号去噪声(分离高频噪声IMF，比如IMF1，然后叠加其它的IMF和残差的大去噪后的信号)； 2、瞬时频率分析(得到IMF后对每个IMF进行希尔伯特变换再提取瞬时频率，也就是所谓的希尔伯特-黄变换，这也是本博文要实践的话题)； 3、从信号的imf分量中进一步提取有用信息（比如：机械故障诊断，通过查看振动信号的imf分量，评估是否有异常；在无人机旋翼参数估测、分类等应用中，基于imf获取旋翼的旋转周期等信息，基于imf构建一些可供分类的特征等等，这也是我本人课题相关）。

如果信号是复数形式的，如何实现EMD？

EMD是处理实数信号的。如果信号是复数，一般有三种处理方法：1、不去管信号的虚部，只拿出信号的实部进行处理即可； 2、处理该复数信号的幅值； 3、分别用EMD分解复数信号的实部和虚部，然后将分解的各个imf结果合并成一个复数imf，然后求其相位，最后得到瞬时频率信息。【我在后文中只给出了当采集的是实数信号情况下的EMD分解和瞬时频率处理结果，这里提及的其它三种情况下的处理只需要在实数信号的情况下做些许代码改动即可，读者可以自行尝试，此外我也在所附代码资料中给出了相关的代码供参考。】

本节的仿真中，我设计了两个chirp信号进行叠加，其参数列表如下：

表3.1 信号参数预设列表

	信号1	信号2
起始频点	20Hz	40Hz
调频斜率	50Hz/s	50Hz/s
采样率	400Hz
采样点数	1024
噪声	SNR = 20dB的高斯白噪声

仿真得到的结果如下：

图3.8 信号时域

如果不对其进行EMD分解，而是直接沿用前面3.1.1的求解过程(直接进行希尔伯特变换得到解析信号，然后求解析信号的相位，再求其频率)，则最终得到的信号频率随时间的变化关系如下：

图3.9 直接进行瞬时频率法求解

其频率并没有出现多个分量，而且其值也与预设的不符合。这种求解方法是不正确的。

对时域信号做EMD分解(Matlab有其自带的emd函数)，我这里预设了三个imf分量，得到的结果如下：

图3.10 EMD分解结果

随后对每个imf分量按照前面单分量的求解流程求其瞬时频率，并叠加到同一幅图中，得到的结果如下：

图3.11 各分量瞬时频率处理结果(原始信号瞬时频率处理结果)

其实效果比较差… 不过也能捕捉到两个主要的频率分量的信息？至于为什么效果比较差，大概是模态混叠的原因。不同的频率成分混入了同一个imf。针对这一问题，有一些改进的方法，比如所谓的集合EMD(eemd),其通过添加高斯白噪声多次分解并平均以抑制模态混叠,具体的原理读者自行阅读相关资料(我也没有去深入研究)，不过Matlab中没有自带的该方法下的函数，需要安装第三方EEMD工具箱才能使用该函数。参考资料[3]中有一个工具箱，我也安装并尝试了其中的eemd函数，但是得到的结果还是不理想：

图3.12 eemd处理后的结果

该函数需要设置噪声强度以及集合次数两个超参，我在仿真中预设为（0.005，100），读者可以尝试其它的参数看看是否有更好的效果。

该工具箱我也一并附在后文的仿真代码链接中了（里面有文档很清晰地告诉你如何安装）。 需要注意的是：因为该工具箱中也包含了emd这个函数，且它工具箱中的emd函数不支持设置输出的imf个数，如果你不想用它这个工具箱中的emd函数，你可以在matlab的搜索路径中将该工具箱中emd函数所在文件移至底端，这样Matlab调用emd函数时才会优先调用其自带的emd函数。