做网站大型,关键词排名点击软件首页,泉州seo网站关键词优,北京seo百度推广引言 本文是信号(瞬时)频率求解与仿真实践专题的第二篇文章#xff0c;在上一篇博文 [1]信号(瞬时)频率求解与仿真实践(1)-CSDN博客中#xff0c;我构建了信号瞬时频率求解的基本框架#xff0c;并且比较详细地讨论了瞬时频率法。这篇博文探讨适用于信号瞬时频率求解的另一种…引言 本文是信号(瞬时)频率求解与仿真实践专题的第二篇文章在上一篇博文 [1]信号(瞬时)频率求解与仿真实践(1)-CSDN博客中我构建了信号瞬时频率求解的基本框架并且比较详细地讨论了瞬时频率法。这篇博文探讨适用于信号瞬时频率求解的另一种方法联合时频分析(对应博文[1]的第3.2节)。 本文将给出几种典型的联合时频分析法对其理论做介绍并实践对仿真得到的回波信号以及雷达的实测数据进行联合时频处理。(各联合时频分析法的具体实现Matlab都有自带的函数本文提供的代码将直接调用这些函数本文重点不是如何自编代码实现这些联合时频算法。) Blog
2025.6.13 本文第一次撰写 目录
引言
目录
一、联合时频分析法概述以及几种典型算法
1.1 概述
1.2 短时傅里叶变换(short time fourier transform, STFT) 及其拓展 1.3 Wigner-Ville分布(Wigner-Ville distribution WVD)及其拓展
1.4 小波变换wavelet Transform, WT
1.5 其它方法
二、本文的仿真信号参数以及实测数据集介绍
2.1 仿真信号参数介绍
2.2 实测数据集数据简介
三、基于各联合时频分析算法的实践
3.1 STFT处理结果
3.2 CVD谱图
3.3 cepstrogram谱图
3.4 WVDSPWVD
3.5 小波变换处理结果
3.6 Capon算法处理结果
四、总结
五、参考资料
六、本文相关仿真代码 一、联合时频分析法概述以及几种典型算法
1.1 概述 联合时频分析法是求解信号瞬时频率的另一种解决思路(还有一种前文介绍的瞬时频率分析法)它不需要求解信号的相位以及做多分量信号的分解而是同时在时间和频率维度处理信号并可以直接给出信号在二维的时间-频率域的能量分布。该方法很早就被用来分析时变频谱其不仅可用于单分量信号也可以直接用于多分量信号的时频分析。 联合时频分析法可以进一步被分为线性时频变换和双线性时频变换。短时傅里叶变换(Short-time Fourier TransformSTFT)就是典型的线性时频变换法常见的双线性时频变换如维格纳-维利分布Wigner-Ville distributionWVD以及在其基础上的各类改进算法。
本章将介绍几种常用的联合时频分析法。 1.2 短时傅里叶变换(short time fourier transform, STFT) 及其拓展 短时傅里叶变换是线性时频率变换的一种也是最常用的联合时频分析法。STFT的数学表达式如下 1-1 式中x(t)是待处理信号w(t-τ)是窗函数。与FFT直接对信号在整个时间维度进行处理相比STFT的本质是通过设置特定的窗函数将信号分割成很多个短的时间段随后对每个时间段进行FFT处理最后将各个时间段处理的结果进行排序从而得到信号的二维时频分布。STFT最大的缺点是需要在时间分辨率和频率分辨率之间做取舍更长的时间段可以获得更好的频率分辨率但是牺牲了时间分辨率反之亦然。在一些特定的应用场景下考虑到当目标速度的变化周期在时间维度发生变化时固定窗函数下可能无法准确捕捉其速度特征也有文章[2]提出了基于信号周期的自适应窗函数下的STFT处理在STFT过程中将信号分成时间间隔不相等的几段进行FFT处理。 虽然STFT下的分辨率受限但是该方法计算简单也不会引入干扰项所以依然是时频分析里面应用最广泛的方法。 我们把STFT处理得到的数据矩阵取幅值后称之为Spectrogram。在STFT得到Spectrogram谱图的基础上还有一些更进一步的处理方法得到新的谱图常用的如Candence Velocity DiagramCVD[3]和Cepstrogram[4],这两谱图都是在STFT的结果上做进一步的处理得到的其关系如下图所示[5] 图1.1 Spectrogram、CVD、Cepstrogram三种谱图的关系 CVD谱图是在Sepctrogram的基础上取对数(自然对数和以10为底的对数都可)然后对每个频率点求时间维度的FFT得到。CVD谱图中包含了时频域中曲线的重复情况提供了关于目标运动部件在时频图中曲线的形状、大小和频率的信息。 Cepstrogram也是在Sepctrogram的基础上取对数(自然对数和以10为底的对数都可)不过随后是对每个时间点求频率维度的IFFT得到Cepstrogram谱图中包含了时频谱的周期性信息。 如上图所示这三种谱图在很多的文献中被用于做微动特征的提取(Feature extraction)或者直接作为样本输入到CNN等深度学习网络中进行训练和目标识别。 Matlab中有自带的函数spectrogram 可以实现STFT功能这里不做过多介绍。 1.3 Wigner-Ville分布(Wigner-Ville distribution WVD)及其拓展 为了提高时频分辨率一些双线性时频变换方法被提出。WVD就是其中最简单以及最有代表性的方法WVD是通过计算信号自相关函数的傅里叶变换得到表达式如下 (1-2) 式中τ是时移变量(在t的基础上做时移)s*表示复共轭是信号在时间t附近的局部自相关其衡量了信号在时间t附近的相似度在此基础上进行傅里叶变换可以得到信号在时间t附近的频率分布。这种方法的时频分辨率远高于线性时频变换但是计算复杂度高且受到交叉项干扰如果信号中有多个频率信息此时对多个频率信号和的WVD并不等于对各个频率信号WVD的和也就是说如果信号中包含多个频点那么WVD处理后会有交叉项而且交叉项对应的幅值要高于各个频率分量如果不做抑制会产生错误的时频谱图影响后续处理。 为了抑制该交叉项的干扰需要进行滤波处理。不同的滤波方法(基于核函数)就有了不同的双线性时频分析方法如下表所示[6] 图1.2 各双线性时频分析法及其核函数 以比较常见的平滑伪wigner-ville分布(Smooth Pseudo WVD,SPWVD)为例其在WVD的基础上引入频域平滑窗h(τ)和时域平滑窗g(s-t)来减少交叉项的干扰 (1-3) 窗函数的引入虽然可以抑制交叉项的干扰但是会影响时间分辨率窗越宽交叉项越少但是时间分辨率也会越低。 1.4 小波变换wavelet Transform, WT 小波变换也是一种常见的时频分析方法不同于傅里叶变换中使用正弦函数作为基函数对信号进行分解小波变换使用可变尺度的基函数(小波)在时频域上对信号进行分析其可以很好地解决STFT无法同时兼顾时间和频率分辨率的缺点。有关小波变换更形象的理解可以参考资料[7]。 小波变换使用母小波ψ(t)通过平移和缩放生成一组小波基 (1-4) 式中a是尺度参数可以控制频率。B是平移参数可以控制时间位置常见的母小波如morlet波、haar波等。连续小波变换公式如下 1-5 小波变换没有分辨率的限制适合精细分析但是该方法计算量大冗余信息多而且选定了母小波后在整个分析过程中将无法更换即使该母小波在全局可能是最佳的但在某些局部可能并不是所以小波分析的基函数缺乏适应性。 1.5 其它方法 A. Capon变换或称为最小方差无失真响应Minimum Variance Distortionless ResponseMVDR Capon变换是一种典型的超分辨率波束形成算法(我在很早前的车载毫米波雷达DOA估计研究[8]中讨论过该算法)将之引入时频分析其实原理与测角是类似的把所采集的时域信号理解成测角时的各通道采集的信号即可。其实现过程可以分为四步
1、信号分段(类似STFT中设置时间窗把原始信号分成很多的小段)
2、对每段信号计算其自相关矩阵R(信号的协方差矩阵)
3、构造频率导向矢量 4、对于每个频点计算该特定时间段和特定频点下的Capon值 通过对每段信号进行前述步骤2到4的处理便可以得到Capon变换下的时频图。Capon方法下相较于STFT、WVD等方法有更高的频率分辨率(时间分辨率没有太多改善)且无交叉项干扰但是该方法需要耗费较大的计算资源而且基于单段信号计算得到的协方差矩阵可能误差较大、求逆不稳定算法的抗噪以及稳定性较差。 B. 联合时频分析算法当然远不止上面几种不过上面几种是我在课题相关的专业文献中遇到使用比较多的几种。 后续如果碰到有其它比较常见的算法我再做补充。 二、本文的仿真信号参数以及实测数据集介绍 本章对待处理的信号做介绍。本文的实践中我预设了一份比较理想的仿真信号、以及使用了一份实测数据集中的信号分别实践各联合时频处理算法。
2.1 仿真信号参数介绍 仿真的信号设置为两个chirp信号的叠加其参数列表如下
表2.1 仿真信号参数列表 参数 Chirp1 Chirp2 起始频率 0 200Hz 调频斜率 100Hz/s -50Hz/s 采样率(双路采样) 600Hz 采样点数 1024 采样时间 1.7067s 计算得到终止频率 170.67Hz 114.67Hz 噪声 SNR20dB的高斯白噪声 这是一个典型的多分量且频率时变信号在本系列上一篇文章[1]中我们使用希尔伯特-黄变换处理过类似的信号但是效果不佳。 信号的时域以及直接FFT处理的结果如下 图2.1 仿真信号时域 图2.2 仿真信号直接进行FFT处理的结果 从上图可以看到信号频段主要位于0至200Hz之间而且在大约110Hz至170Hz频段内能量比较大这刚好是我们设置的两信号分量的频段重叠区域。综上可知信号的频段表现与预设信号参数是吻合的。 2.2 实测数据集数据简介 实测数据我选取了自己课题方向相关的一份雷达探测无人机的数据集。(我在之前的博文[9]中做过雷达测量人体运动的回波时频图处理感兴趣的读者也可以下载博文[9]中的数据集进行处理)。 本文使用的这份数据集读者可以从信号处理学报的官网上找到下载链接[10]数据集与这份数据集相对应的文章是[11],文中有关于该数据集的详细介绍实验对象、实验场景、数据文件使用说明等作者在文中还给出了一些基于本数据集的处理示例。本文就不对该数据集做细节的介绍了本文后续的实践只是取出其中一部分数据做时频处理(在后文的代码链接中我一并附上了后续章节使用的(从其原始数据中解析出来的原始回波数据的)二维数据矩阵.mat格式读者直接load即可当然我也附上了数据解析和预处理相关的代码。 下表为该数据集所使用雷达的主要参数
表2.2 实测数据集雷达参数列表 参数 值 载频 3.13344GHz 收发通道数 1T8R 工作体制 调频连续波 发射带宽 120MHz 脉冲宽度 2.0835微妙 脉冲重复周期 533.3或1000微秒 中频采样频率 491.52MHz 本文后续的实践部分我分别选取了无人机悬停选用其数据集Data3文件夹下的前八个数据文件和无人机运动两种飞行状态下选用其数据集Data1文件夹下的前八个数据文件各一份数据进行时频处理。 其单份小文件下的数据量是8*128*81928对应的是8个接收通道8192对应单个chirp下距离维采样点数128对应这份数据文件下包含了128个chirp。 【虽然其真值文件中给出了目标的角度但是其数据说明文档中没有给出接收阵列的排布方式所以我在做处理时没有使用数字波束形成的方法去增加目标的SNR,而是直接只使用其中一个通道的数据】此外所选用的8个小数据文件其实前后是时间连续的所以这八个数据文件可以在chirp维度连接起来。综上经过数据解析和排布得到了两个(分别对应无人机悬停和飞行)大小为1024*8192 的原始回波数据矩阵(这两份数据我保存为.mat格式一并附在代码链接中读者可以直接load后使用而不用跑数据解析部分的代码) 随后对该矩阵进行匹配滤波(距离维压缩) --- MTI去除静态杂波 ---- 多普勒维FFT处理得到RD图 ---- 找到目标所在距离门 ---- 取出MTI后的数据矩阵中目标所在距离门下的数据进行时频处理 注实践中发现对于无人机悬停的情况(此时机身运动速度为0)即便进行MTI处理也还是可以从RD图中看到无人机回波而不做MTI处理后续时频处理后机身的能量很大会淹没很多的旋翼微动信息。【读者可以基于我提供的代码做更一些对比分析(直接把代码中MTI处理的部分关闭即可)】 1. 无人机悬停状态下所获取数据的预处理结果 图2.3 回波时域幅值(一共1024个脉冲) 图2.4 回波时域(取其中一个脉冲) 随后对其进行匹配滤波处理得到 图2.5 距离压缩后的结果 同样取其中一个脉冲画图 图2.6 距离压缩后的结果(取其中一个脉冲) 从图中可以看到杂波很强并没有明确的目标位置。对距离压缩后的结果做MTI有关静止目标滤除的方法读者可以参考我之前的博文[12],得到结果如下 图2.7 距离压缩后的结果MTI后 此时可以看到在距离索引[1400 1900]之间能量比较高可以认为目标就在该区间内。进一步地对MTI后的结果再做脉冲维的FFT处理得到距离-多普勒图(RD图)如下 图2.8 回波数据矩阵对应的RD图 在这份数据集的说明文档中Data3对应的实验场景为无人机处于悬停状态其与雷达之间的斜距约为212.83m这与上图的结果基本吻合。 【注图中之所以在速度维度有周期性的波峰这是无人机旋翼旋转产生的所谓HERM line当叶片旋转时叶片会对雷达电磁波产生周期性调制本质就是微多普勒信号。如果解调时间窗很长(上面我们用到了全部1024个脉冲)就会产生上图中的HERM line包括很多的谐波如果解调时间窗比较短会出现所谓的blade flash(在后文的时频处理中会看到)我们可以从时频图中看到叶片产生的类似正弦信号的调制波形。 进一步地我们可以从这两种图中获取很多的无人机参数信息如叶片旋转速率、叶片长度、旋翼个数等参数。不过这些内容不在本文的讨论范围之内】 在后文的时频处理中我选取了MTI处理后的距离索引为[1400 1500]区间内的数据做处理。 2. 无人机运动状态下所获取数据的预处理结果 图2.9 回波时域幅值(一共1024个脉冲) 图2.10 回波时域(取其中一个脉冲) 随后对其进行匹配滤波处理得到 图2.11 距离压缩后的结果 图2.12 距离压缩后的结果(取其中一个脉冲) MTI处理后得到 图2.13 距离压缩后的结果MTI后 进一步地得到距离-多普勒图(RD图)如下 图2.14 回波数据矩阵对应的RD图 在这份数据集的说明文档中Data1对应的实验场景为无人机处于由近及远以S型轨迹飞行的状态速度约5.8m/s其与雷达之间的斜距约为250.41m这与上图的结果基本吻合(速度的绝对值有一些差异可能是其真值测量不准确)。 在后文的时频处理中我选取了MTI处理后的距离索引为[1550 1650]区间内的数据做处理。 三、基于各联合时频分析算法的实践 本章给出各联合时频分析算法对第二章所介绍的仿真信号以及实测数据的处理结果。 (相关的实现代码可参考第六章给的代码链接)。 本章的处理结果图较多在实测数据处理部分我后文提供的代码都会同时输出纵轴分别为频率和速度的两幅图不过在本文给的结果中我只选取了其中一幅(而且也没有统一)。 本章的处理结果不一定是最优的各类时频处理方法中有很多的参数可调节读者可以基于后文提供的代码自行设置和优化参数。 3.1 STFT处理结果
A. 仿真数据处理结果 STFT的处理参数为时间窗长64前后时间窗的重叠率为90%加汉宁窗。 图3.1 STFT对仿真数据的处理结果 结果是符合2.1节中的预设参数的。 B实测数据处理结果(目标悬停) 预设时间窗长为16前后窗的重叠率为0.9得到结果如下 图3.2 STFT对实测数据(目标悬停)的处理结果blade flash图 还是可以看到一些由叶片旋转引起的周期性的变化此外通过把窗长增加到128可以看到HERM line图 图3.3 STFT对实测数据(目标悬停)的处理结果HERM line图 这些沿速度维排布的能量较大的线条就是所谓的HERM line图这与我们在前面得到RD-mat类似。 C实测数据处理结果(目标运动) 预设时间窗长为16前后窗的重叠率为0.9得到结果如下 图3.4 STFT对实测数据(目标悬停)的处理结果blade flash图 可以看到速度为-5到-10之间能量很强对应机身的运动。在-15周围也有一些周期性出现的能量强点这些对应的就是旋翼。通过把窗长增加到128可以看到HERM line图 图3.5 STFT对实测数据(目标悬停)的处理结果HERM line图 3.2 CVD谱图
A. 仿真数据处理结果 经由图3.1的结果处理而来 图3.6 对仿真数据处理得到的CVD图 似乎看不到什么有用的信息。 B实测数据处理结果(目标悬停) 设置STFT的窗长为16在STFT的结果上处理得到(这里采用取自然对数的方式实现) 图3.7 对目标悬停状态下回波处理得到的CVD图 结果暂不做评价。 读者可以基于后文提供的代码 尝试设置不同的STFT参数以及取以10为底的对数或者对整个STFT的结果先做归一化处理再进行1.2所介绍的CVD处理。看看是否有其它更好的结果。 C实测数据处理结果(目标运动) 图3.8 对目标运动状态下回波处理得到的CVD图 3.3 cepstrogram谱图
A. 仿真数据处理结果 经由图3.1的结果处理得到 图3.9 对仿真数据处理得到的Cepstrogram图 似乎看不到什么有用的信息。 B实测数据处理结果(目标悬停) 设置STFT的窗长为16在STFT的结果上处理得到(这里采用取自然对数的方式实现) 图3.10 对目标悬停状态下回波处理得到的Cepstrogram图 结果暂不做评价。 C实测数据处理结果(目标运动) 设置STFT处理时的窗长为128得到结果如下 图3.11 对目标运动状态下回波处理得到的Cepstrogram图 3.4 WVDSPWVD
A. 仿真数据处理结果 图3.12 对仿真数据做WVD处理的结果 虽然可以看到预设频率下的信号但是整个时频图中有比较强烈的干扰。随后预设时间和频率窗长都为129得到SPWVD的处理结果如下 图3.13 对仿真数据做SPWVD处理的结果 可以看到交叉干扰项确实得到了有效抑制。读者可以基于给的代码设置其它长度的窗长看看效果。 B实测数据处理结果(目标悬停) 图3.14 对目标悬停状态下的回波数据做WVD处理的结果 WVD无法输出负频带的信息这部分信息会被耦合到正频带中。 设置时间和频率窗的长度都是127得到SPWVD的处理结果如下 图3.15 对目标悬停状态下的回波数据做SPWVD处理的结果 C实测数据处理结果(目标运动) 图3.16 对目标运动状态下的回波数据做WVD处理的结果 设置时间和频率窗的长度都是127得到SPWVD的处理结果如下 图3.17 对目标运动状态下的回波数据做SPWVD处理的结果 这两种方式处理得到的时频图比较接近STFT处理的结果。 3.5 小波变换处理结果 本节的处理中都使用’morse’小波基函数。
A. 仿真数据处理结果 图3.18 CWT变换对仿真数据的处理结果 可以看到与STFT变换类似的效果但是动态范围低很多。(注意频率范围的变化) B实测数据处理结果(目标悬停) 图3.19 CWT变换对实测数据(目标悬停)的处理结果 C实测数据处理结果(目标运动) 图3.20 CWT变换对实测数据(目标运动)的处理结果 还是可以看到在速度为5m/s左右有一条能量较大的曲线。 3.6 Capon算法处理结果 本节的仿真中都预设窗长为64前后时间段移动16点频率点数为256点。 A. 仿真数据处理结果 图3.21 Capon算法对仿真数据的处理结果 效果与STFT的处理结果类似不过SNR要低一些。 B实测数据处理结果(目标悬停) 图3.22 Capon算法对实测数据(目标悬停)的处理结果 C实测数据处理结果(目标运动) 图3.23 Capon算法对实测数据(目标运动)的处理结果 可以看到其准备捕捉到了目标的运动速度不过没有旋翼转动引起的速度调制等细节。 四、总结 时频图的获取往往并不是我们的最终目的时频图是一种很好的、可用于分析雷达目标的媒介。一般来说获取时频图主要有三个目的一是从图中提取所需的目标参数信息(比如无人机旋翼数量、叶片长度、叶片旋转频率等参数)二是从时频图中提取一些特征以供给后续的传统机器学习算法三是将时频图直接给到深度学习网络中进行训练并完成目标识别等任务。 本文作为信号(瞬时)频率求解与仿真实践话题下的第二篇文章主要介绍并实践了几种典型的联合时频处理算法。 五、参考资料
[1] 信号(瞬时)频率求解与仿真实践(1)-CSDN博客
[2] D. B. Herr and D. Tahmoush, Data-Driven STFT for UAV Micro-Doppler Signature Analysis, 2020 IEEE International Radar Conference (RADAR), Washington, DC, USA, 2020, pp. 1023-1028, doi: 10.1109/RADAR42522.2020.9114726.
[3] S. Björklund, T. Johansson, and H. Petersson, “Evaluation of a micro Doppler classification method on mm-wave data,” in Proc. IEEE Radar Conf., May 2012, pp. 934–939.
[4] A. M. Noll and M. R. Schroeder, “Short-time ‘cepstrum’ pitch detection,” J. Acoust. Soc. Amer., vol. 36, no. 5, p. 1030, 1964, doi: 10.1121/1.2143271.
[5] A. Hanif, M. Muaz, A. Hasan and M. Adeel, Micro-Doppler Based Target Recognition With Radars: A Review, in IEEE Sensors Journal, vol. 22, no. 4, pp. 2948-2961, 15 Feb.15, 2022, doi: 10.1109/JSEN.2022.3141213.
[6] V. C. Chen, F. Li, S. . -S. Ho and H. Wechsler, Micro-Doppler effect in radar: phenomenon, model, and simulation study, in IEEE Transactions on Aerospace and Electronic Systems, vol. 42, no. 1, pp. 2-21, Jan. 2006, doi: 10.1109/TAES.2006.1603402.
[7] https://zhuanlan.zhihu.com/p/22450818
[8] 车载毫米波雷达DOA估计综述-CSDN博客
[9] 雷达一维成像基于数据集的实践-CSDN博客
[10] 数据集
[11] 林钱强, 秦正阳, 蒋李兵, 等. 简单野外背景条件下地面雷达对“低慢小”无人机探测数据集[J]. 信号处理, 2024, 40(11): 2095-2104. DOI: 10.12466/xhcl.2024.11.014.
[12] 毫米波雷达信号处理中的静止目标(静态杂波)滤除问题-CSDN博客 六、本文相关仿真代码 本文相关的代码、可以直接load的实测数据原始数据集很大感兴趣的读者可以自行前往信号处理学报[10]中下载等资料一并打包在如下链接中读者可自行下载。
https://download.csdn.net/download/xhblair/90999073
信号(瞬时)频率求解与仿真实践(2)博文相应的资料资源-CSDN文库
