现代数字信号处理课程报告

问题：需要利用4种时频分析方法对给定的6个信号文件进行处理，给出每种分析方法的结果图。由于时频分析方法众多，选取的用于MATLAB编程的4种时频分析方法分别为：短时傅立叶变换（STFT）、连续小波变换、Wigner-Ville分布、希尔伯特黄变换（HHT）。短时傅立叶变换：傅立叶分析一直被认为是最完美的数学理论和最实用的方法之一。但是用傅立叶分析只能获得信号的整个频谱，而难以获得信号的局部特性，特别是对于突变信号和非平稳信号难以获得希望的结果。为了克服经典傅立叶分析本身的弱点，1946年Gabor提出了短时傅立叶变换，其基本思想是首先利用窗函数截取信号，假设信号在各个窗内都是平稳的，再对第一个窗内的信号采用傅立叶变换进行分析，以确定此时窗内存在的频率成分，然后沿着信号时间方向移动窗函数，得到频率随时间的变换关系，即所需要的时频分布。解题思路与流程：分别读取6个信号，对读取的信号，计算其长度，从而选取窗口长度为100，因此重叠宽度为50，从而计算需要分段的数目，进而计算做fft的点数，此处采取的窗函数为矩形窗，沿着信号时间方向移动窗函数，从而求出每段的fft，拼接后画出时频分析二维与三维图。6个信号的图如图1-6所示。a.NACB.BHR三维b.NACB.BHR二维图1NACB.BHRa.NACB.BHT三维b.NACB.BHT二维图2NACB.BHTnfNACB.BHR信号STFT200400600800100012001400-0.4-0.3-0.2-0.100.10.20.30.40.5nfNACB.BHT信号STFT5005506006507007508008509009501000-0.8-0.6-0.4-0.200.20.40.60.8a.NACB.BHZ三维b.NACB.BHZ二维图3NACB.BHZb.YHNB.BHR三维b.YHNB.BHR二维图4YHNB.BHRc.YHNB.BHT三维b.YHNB.BHT二维图5YHNB.BHTd.YHNB.BHZ三维b.YHNB.BHZ二维图6YHNB.BHZnfNACB.BHZ信号STFT200400600800100012001400-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5nfYHNB.BHR信号STFT200400600800100012001400-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5nfYHNB.BHT信号STFT200400600800100012001400-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5nfYHNB.BHZ信号STFT200400600800100012001400-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5连续小波变换：由于STFT采用的滑动窗函数一经选定就固定不变，故决定了其时频分辨率固定不变，不具备自适应能力，因此短时傅立叶变换没有从根本上解决傅立叶分析的固有问题。小波变换的诞生，正是为了克服经典傅立叶分析本身的不足。小波变换的基本思想是通过伸缩平移运算对信号逐步进行多尺度细化，最终达到高频处时间细分和低频处频率细分，从而聚焦到信号的任意细节。解题思路与流程：由于小波分析方法是一种窗口面积固定但其形状可以改变的方法，因此我们弹性的时频窗，将信号分为了128个尺度，不同的尺度对应不同的频率，频率高的地方，需要时间间隔小，因此高频区会有一个较高的时间分辨率，而低频的地方会有较高的频率分辨率，这样就保证了时间分辨率和频率分辨率在各自需要的范围都达到很高的精度。在本次实验中选取总尺度为128，选取sym2函数为小波函数，计算出小波基的中心频率，根据中心频率与总尺度计算出各个尺度划分，从而求出小波系数，进而画出如图7-12的小波时频图。a.NACB.BHR三维b.NACB.BHR二维图7NACB.BHRa.NACB.BHT三维b.NACB.BHT二维图8NACB.BHT时间t/s频率f/HzNACB.BHR信号sym2小波时频图20040060080010001200140016000.050.10.150.20.250.30.350.40.450.5123456789x10-6时间t/s频率f/HzNACB.BHT信号sym2小波时频图20040060080010001200140016000.050.10.150.20.250.30.350.40.450.524681012x10-6a.NACB.BHZ三维b.NACB.BHZ二维图9NACB.BHZa.YHNB.BHR三维b.YHNB.BHR二维图10YHNB.BHRa.YHNB.BHT三维b.YHNB.BHT二维图11YHNB.BHTa.YHNB.BHZ三维b.YHNB.BHZ二维图12YHNB.BHZ时间t/s频率f/HzNACB.BHZ信号sym2小波时频图20040060080010001200140016000.050.10.150.20.250.30.350.40.450.524681012141618x10-6时间t/s频率f/HzYHNB.BHR信号sym2小波时频图20040060080010001200140016000.050.10.150.20.250.30.350.40.450.50.20.40.60.811.21.41.61.822.2x10-5时间t/s频率f/HzYHNB.BHT信号sym2小波时频图20040060080010001200140016000.050.10.150.20.250.30.350.40.450.50.20.40.60.811.21.41.61.82x10-5时间t/s频率f/HzYHNB.BHZ信号sym2小波时频图20040060080010001200140016000.050.10.150.20.250.30.350.40.450.50.511.522.5x10-5Wigner-Ville分布：Wigner-Villedistribution(WVD)方法，是由Wigner1932提出并由Ville于1948年引入信号分析领域。WVD可看作信号能量在联合的时间和频率域中的分布，其功能是分析时间序列中的频率成分及其能量，是分析非平稳和时变信号的重要工具。WVD基本思想为信号中心协方差函数的傅立叶变换，它具有许多优良的性能，如对称性、时移性、组合性、复共轭关系等，不会损失信号的幅值与相位信息，对瞬时频率和群延时有清晰的概念。其不足是不能保证非负性，尤其是对多分量信号或具有复杂调制规律的信号会产生严重的交叉项干扰，这是二次型时频分布的固有结果，大量的交叉项会淹没或严重干扰信号的自项，模糊信号的原始特征。解题思路与流程：对于x(t)为一连续时间信号，则公式（1）称为信号x(t)的自WVD。𝑊𝑉𝐷𝑥(𝑡,𝜔)=∫𝑥(𝑡+𝜏/2)𝑥∗(𝑡−𝜏/2)exp⁡(−𝑗𝜔𝜏)𝑑𝜏（1）式（1）中，𝑥∗(𝑡)为𝑥(𝑡)的复共轭，𝜏为频率轴划分区间，值为2的n次方，因此对x(t)信号进行处理时，可以先计算信号的𝑥(𝑡+𝜏/2)𝑥∗(𝑡−𝜏/2)，然后对结果取fft，从而画出如图13-18所示时频分析图。图13NACB.BHR图14NACB.BHT图15NACB.BHZ图16YHNB.BHR图17YHNB.BHT图18YHNB.BHZHHT变换：HHT就是Hilbert-HuangTransform（希尔伯特黄变换），是由美国航天航空局黄鳄教授于1998年提出。HHT方法从信号自身特征出发，用经验模态分解（EMD）方法把信号分解成一系列的本征模态函数（IMF），然后对这些IMF分量进行Hilbert变换，从而得到时频平面上能量分布的Hilbert谱图，打破了测不准原理的限制，可以准确地表达信号在时频面上的各类信息。HHT的功能是可以分析并提取时间序列中的主要频率（即经验模态）成分。其特点在于可以提取（而不仅仅是分析）信号中的主要成分。HHT从数学理论上很完美，但在应用时常常受到攻击：其所谓经验模态缺乏物理本质或意义。解题思路与流程：根据信号的长度，先对信号进行经验模态分解（Empiricalmodedecomposition,EMD），此时是利用时频工具箱中的emd函数，得出本征模态函数（intrinsicmodefunction,IMF），再对本征模态函数进行希尔伯特变换，从而过进一步得该信号的希尔伯特谱、时频能量谱等，以便对信号进行分析，对非线性及非平稳信号有较好的分析和处理效果。流程如下：X(t)本征模态函数（IMF）时频能量谱EMDHilbert变换利用MATLAB编程实现，得到如图19-24所示的三维时频分析图。图19NACB.BHR图20NACB.BHT图21NACB.BHZ图22YHNB.BHR图23YHNB.BHT图24YHNB.BHZ结果分析：根据图1.b-6.b（利用STFT处理6个信号二维时频图）可知，图1.b信号的时间分辨率较高，但是频率分辨率能力差；图2.b信号的频率分辨率较高，而时间分辨率较差，这是因为图1.b在利用短时傅立叶变换时，选取的窗口长度（Nw=20）过短，导致有较高的时间分辨率，但是频率分辨率能力差；图2.b选取的窗口长度（Nw=1000）过长，导致有较高的频率分辨率，但是时间分辨率能力差。根据图7.b-12.b（利用连续小波变换处理6个信号二维时频图）可知，连续小波变换能给出比较好的时间精度，在低频区也有很好的频率精度，但是在高频区频率分辨能力较弱。解决了STFT时间分辨率差的问题。根据图13-18（利用WVD处理6个信号二维时频图）可知，WVD方法在低频区有很好的时间分辨率和频率分辨率，即WVD有很好的时频聚焦性。但是式（1）中，𝑥(𝑡+𝜏/2)𝑥∗(𝑡−𝜏/2)导致信号x(t)出现了两次，而且WVD里不含有任何窗函数，因此必然会产生交叉项，而这些交叉项会在信号的各个分量之间产生虚假信号，导致分析不准确。根据图19-24（利用HHT处理6个信号二维时频图）可知，HHT变换有很高的时间分辨率和频率分辨率，然而EMD方法还不成熟，不知道为什么求出来的二维时频谱只有离散的一些点（只画出了其中一个信号，如下图）。猜想是因为在低频区出现了一些不合理的频率成分，掩盖了低能量的频率成分。程序代码：STFT：clearall；clc;closeall;fid=fopen('NACB.BHR','r','ieee-le');A=fread(fid,[70,1],'float32');B=fread(fid,[40,1],'int32');C=char(fread(fid,[1,192],'char'));x1=fread(fid,'float32');A(A==-12345.0)=NaN;B(B==-12345)=NaN;fclose(fid);figure(1);timefreq(x1,20,'rec');title('NACB.BHR信号STFT');fid=fopen('NACB.BHT','r','ieee-le');A=fread(fid,[70,1],'float32');B=fread(fid,[40,1],'int32');C=char(fread(fid,[1,192],'char'));x2=fread(fid,'float32');n频率f/Hz二维hilbert时-频-幅度谱200400600800100012001400160000.050.10.150.20.250.30.350.40.450.50123456x10-6A(A==-12345.0)=NaN;B(B==-12345)=NaN;fclose(fid);figure(2);timefreq(x2,1000,'rec');title('NACB.BHT信号STFT');fid=fopen('NACB.BHZ','r','ieee-le');