matlab画时域和频谱图_信号的频域分析及matlab实现

随机振动信号分析方法总结

信号处理（信号滤波、时频域分析、神经网络、寿命预测）

一、时域分析

时域分析特征包括均值、方差、峭度、峰峰值等；

振动信号降噪结果分析： 对于去噪效果好坏的评价，常用信号的信噪比（SNR）、估计信号同原信号的均方根误差（RMSE）来判断。SNR 越高则说明混在信号里的噪声越小，否则相反。RMSE的计算值越小则表示去噪效果越好。 信噪比定义：

均方根误差定义：

二、频域分析

三、时频联合域分析（Joint Time-Frequency Analysis，JTFA）

即时频分析，通过对时变非平稳信号提供时域与频域的联合分布信息，从而清楚地描述出时间和频率的相互变化关系。在传统的信号处理上，常使用傅里叶变换及其反变换进行时频分析。但是其作为一种信号的整体变换，不具备时间和频率的定位功能，更适用于平稳信号的分析。而现实中采集的信号多为非平稳信号，由于其频率随时间变化较大，因此分析方法必须能够准确地反映出信号的局部时变频率特性，若使用传统的方法则很难对信号进行更好的分析，因此需要把整体谱推广到局部谱中。时频分析方法是将一维时域信号映射到二维时频平面，全面地反映信号的时频联合特征，通过设计时间和频率的联合函数，从而描述信号在不同时间和频率的能量密度与强度。目前常用的时频分析方法有小波变换、短时傅里叶变换、Wigner-Ville 分布、Hilbert-Huang 变换等。

1. HHT(Hulbert-Huang变换)

HHT 分析法由经验模态分解（Empirical Mode Decomposition，EMD）和 Hilbert 谱分析（Hilbert Spectrum Analysis，HAS）两部分构成。其处理非平稳信号的基本过程如下：首先使用 EMD 将采集到的各类信号分解为满足相应条件的若干固有模态函数（Intrinsic Mode Function，IMF)；然后对每一阶的 IMF 分量进行 Hilbert 变换，得到相应的 Hilbert 谱；最后将这些分量的 Hilbert 谱进行汇总，即可得到原信号的 Hilbert 谱。

1.1 图解经验模态分解(EMD)

这部分转自简书：一位学有余力的同学

1.1.1 经验模态分解是做什么的

**经验模态分解(Empirical Mode Decomposition，EMD)**是由美国工程师黄锷于1998年提出的一种信号的时频分析方法，这里的信号指的是时序信号。

常见的时序信号处理方法可以分为三类：时域、频域和时频域。时域分析特征包括均值、方差、峭度、峰峰值等；频域特征包括频率、能量等；而时频域分析有小波变换等。经验模态分解就属于一种时频分析方法。

黄锷认为所有的信号都是由有限个**本征模函数(Intrinsic Mode Function, IMF)**组成。IMF分量包含了原信号的不同时间尺度的局部特征信号。经验模态分解法能使非平稳数据进行平稳化处理，然后进行希尔伯特变换获得时频谱图，得到有物理意义的频率。

这和快速傅里叶变换(Fast Fourier Transform, FFT)有些像，FFT假设所有信号都是由很多周期性的正弦信号组成，这些信号有着不同的幅频和相位。使用FFT可以将时域信号转换到频域，但EMD分解后的信号还在时域，并且它没有假设信号是周期的且由很多基本的正弦信号组成。

1.1.2 经验模态分解的使用条件

但是EMD的使用存在一些限制条件： ⑴函数在整个时间范围内，局部极值点和过零点的数目必须相等，或最多相差一个； ⑵在任意时刻点，局部最大值的包络（上包络线）和局部最小值的包络（下包络线）平均必须为零。

第一条什么意思呢，看看下面的图就明白了，它只能是下面这种情况：

这里面局部极值点有三个，而过零点有四个，相差一个是符合条件的。而不是这种情况：

第二条的意思就是这个信号的包络线必须是上下对称的，即信号的均值为0。

1.1.3 EMD分解步骤

假如我们有如下信号，它是由频率为1hz和4hz的正弦信号叠加而成：

叠加之后的信号:

第一步：找到信号的极大值与极小值极

第二步：使用三次样条插值法拟合极大值与极小值的包络线

第三步：将两条极值曲线平均获得平均包络线

第四步：用原始信号减去均值包络线

这样我们就得到了第一个IMF，是不是和4hz的信号很像，但是和真实的4hz信号还有一些误差，比如信号的首尾两端的幅值突然增加，这是由于三次样条差值所产生的误差。我们发现得到的这个IMF同样满足EMD的两个条件，我们可以对该IMF作为输入信号从第一步开始计算第二个IMF，直到最终得到的信号是一个常数、单调或者只有一个极值为止。

1.2 EMD分解存在问题及优化

1.2.1 固有模态分解

在实际信号采集的过程中，大多数信号都是复杂信号，含有多个频率成分且在任意时刻的数据可能包含多个震荡模式，难以准确分析测量，因此考虑将信号分解成一系列的单频率分量信号，即把数据分解成固有模态函数（IMF）。 Huang 认为任何信号都可由若干个固有模态函数相互叠加得到，即任何信号都可分解成一系列的局部均值为零的对称信号。

也就是说，一个函数进行 EMD 分解后所获得的每个 IMF 必须具有以下两个条件，才能具有实际的意义： ① 在所有时间内，函数的局部极值点个数和必须与过零点数差值在[-1，1]范围内。 ② 在任一时刻点，由极大值和极小值组成的上下包络线的平均值必须为零。满足上述两个条件之后得到的 IMF 分量能清晰地反映出信号内部的固有频率成分，因此理论上每阶 IMF 分量在各周期内只含一种频率成分，且不会出现模态混叠现象。

1.2.2 模态混叠现象分析

在进行 EMD 分解时，首先需要通过查找极值点来构造包络曲线。信号中存在异常事件时，会对极值点的选取产生很大的影响，使极值点分布不均匀，从而导致求取的包络线为异常事件的包络与正常信号包络的组合。经过该包络计算出的均值而获得的 IMF分量也就随之包含了信号的固有模式和异常事件，或者包含了相邻特征时间尺度的固有模式，从而产生了模态混叠现象，造成 IMF 分量不精确。

简单来说，模态混叠可以理解为一个 IMF 分量中包含差异极大的特征时间尺度或者相近的特征时间尺度被分布在不同的 IMF 分量中，从而导致相邻的两个 IMF 波形混叠，难以辨认的现象。模态混叠产生的主要原因有： ① 当振动信号中混有间断信号，即高频小幅值的信号时，会导致 EMD 分解出现模态混叠现象。 ② 信号中含有大量的噪音信号。 ③ 在多个频率共存的混合信号中，如果不能将某几个频率信号正确分离，也会产生模态混叠现象。

1.2.3 模态混叠抑制方法

分析可知，模态混叠问题的主要原因是在求取包络过程中局部极值在很短的时间间隔内发生多次跳变。针对 EMD 分解的不足法国的 Handrin 等人用 EMD 对白噪声分解后的结果进行统计，提出了一种基于辅助噪声数据分析的改进的 EMD 方法，即集合经验模态分解法（Ensemble Empirical Mode Decomposition，EEMD）【详见Mr.括号大佬—类EMD的“信号分解方法”及MATLAB实现（第一篇）——EEMD】。在进行实验时，利用白噪声频谱均匀分布的特性，在待分析信号中加入白噪声，这样不同时间尺度的信号可以自动分离到与其相适应的参考尺度上去。 EEMD 是一种通过添加噪声进行辅助分析的方法。其原理是给原信号中首先加入白噪声，当附加的白噪声均匀分布在整个时频空间时，该时频空间就由滤波器组分割成的不同尺度成分组成。当信号加上均匀分布的白噪声背景时，不同尺度的信号区域将自动映射到与背景白噪声相关的适当尺度上去，以此来补充一些缺失的尺度，在信号分解中具有良好的表现。由于零均值噪声的特性，噪音经过多次平均计算后会相互抵消，这样最终集成平均的计算结果就能直接视作最终结果。EEMD 分解的流程图：

EEMD 分解过程的主要步骤如下： ① 对于采集到的信号，首先加入具有正太分布的白噪声； ② 将加入白噪声后的混合信号作为一个整体，然后进行 EMD 经验模态分解，得到各 IMF 分量； ③ 重复步骤①和②，每次加入不同的白噪声序列； ④ 将最终得到的多组 IMF 分量进行集成平均处理，计算结果即为最终得到的经优化过的 IMF 分量。

2. 边际谱与傅里叶谱的比较

Matlab论坛cwjy 意义不同：边际谱从统计意义上表征了整组数据每个频率点的累积幅值分布，而傅里叶频谱的某一点频率上的幅值表示在整个信号里有一个含有此频率的三角函数组分。 作用不同：边际谱可以处理非平稳信号，如果信号中存在某一频率的能量出现，就表示一定有该频率的振动波出现，也就是说，边际谱能比较准确地反映信号的实际频率成分。而傅里叶变换只能处理平稳信号。
CSDN括号先森 在傅里叶谱中，在某一频率上存在着能量意味着具有该频率的正弦或余弦波存在于信号的整个持续时间内；而在边际谱中，在某一频率上存在着能量意味着具有该频率的波在信号的整个持续时间内某一时刻出现的可能性较高。因此在一定程度上，Hilbert边际谱具有一定的概率意义，Hilbert边际谱可以看作是一种加权的联合 幅值-频率-时间分布，而赋予每个时间-频率单元的权重即为局部幅值，从而在Hilbert边际谱中，在某一频率上存在着能量就意味着具有该频率的振动存在的可能性，而该振动出现的具体时刻在Hilbert谱中给出。

简单来说，傅里叶谱和边际谱有一定的相关性，但是在处理非平稳信号时，更适合使用边际谱，因为“傅里叶变换为了在数学上拟合原始数据的非平稳波形，不得不引入大量高频的’伪’谐波分量，这会导致傅里叶谱对低频能量的低估

边际谱程序：

代码语言：javascript复制

bjp2 = sum(hs,2) * 1/fs;

代码语言：javascript复制

%网上的程序和  bjp2 = sum(hs,2) * 1/fs;结果一样
for k=1:size(hs,1)
    bjp(k)=sum(hs(k,:)*1/fs);
end
plot(f,bjp)

边际谱和采样频率，数据量有关，此处截取相同数据量；如果数据量时长、采样率不同时：

Matlab论坛宋老师 1，不同的采样率得到边际谱的频率范围不一样，可以用以下两种方法进行调整： (1)把得到的原始信号都通过resample函数调整到相同的采样频率； (2)把边际谱的横坐标不用实际频率，而用归一化频率。 2，不同时长得到的边际谱幅值不同，则是否可以在得到边际谱后除以信号长度N。

全部程序：

代码语言：javascript复制

fs=100;
data = xlsread('C:UsersAdministratorDesktopleosun振动paperZHENDONG数据不同冰型6000dates.xlsx',15);
x=data(:,6);
N = length(x);
L=N;
%% EMD和HT
figure()
emd(x);
[imf,residual,info]=emd(x,'Interpolation','pchip','Display',0);
% figure()
% hht(imf,fs);
% 横轴表示时间、纵轴表示频率，颜色表示能量
[hs, f, t, imfinsf, imfinse] = hht(imf,fs);
% hs——信号的希尔伯特谱（Hilbert Spectrum ）
% f——信号的频率向量（Frequency vector of signal）
% t——信号的时间向量（Time vector of signal）
% imfinsf——每个imf的瞬时频率（instantaneous frequency of each imf）
% imfinse——每个imf的瞬时能量（instantaneous energy of each imf）
[m,n] = size(hs);
l = size(imf,2);

% 画出所有imf函数。
% figure()
% for i = 1:l
%     plot(1:n,imfinsf(:,i));
%     hold on;
% end

%% 由于HHT变换，得到的时频谱hs是错误排列，所以需要重建时频谱
% 重建思路：HHT变换得到每一条imf的瞬时频率向量，和对应的瞬时能量向量
% f = (0:m-1)/(m-1)*(fs/2)，这个是记录每一个点的记录频率，反向操作，就可以求出瞬时频率对应的k(0<=k<=100)
k = imfinsf*(m-1)*2/fs;
k = round(k);

% f_k是一个记录频率次序的矩阵，当次序k<0时，f_k=0，对应的能量为0;
% 次序>=0的，能量为原来的值，
f_k = k;
f_k(find(f_k<0)) = 0;

k_se = zeros(size(k));
k_se(find(k>=0))=1;
se = imfinse .* k_se; % 得到转换后的，能量矩阵

% 得到时频谱
spp = zeros(m,n);

for i = 1:n
    spp(f_k(i,:) 1,i) = spp(f_k(i,:) 1,i)   se(i,:)';
end




% 得到边际谱
spp_fd = sqrt(spp);

figure;

Y = fft(x);
P2 = abs(Y/L);
P1 = P2(1:L/2 1);
P1(2:end-1) = 2*P1(2:end-1);

ff = fs*(0:(L/2))/L;

subplot 311;plot(ff,P1,'linewidth',1.5);xlabel('频率/(Hz)','fontsize',12);ylabel('幅值/（m·s-2）','fontsize',12);title('N1平整冰垂向-频谱','fontsize',12);



subplot 312;
bjp2 = sum(hs,2) * 1/fs;
plot(f,bjp2);
xlabel('频率/(Hz)','fontsize',12);ylabel('幅值/（m·s-2）','fontsize',12);
title('垂向-边际谱','fontsize',12);


subplot 313
bjp1 = sum(spp_fd,2) * 1/fs;
plot(f,bjp1);
xlabel('频率/(Hz)','fontsize',12);ylabel('幅值/（m·s-2）','fontsize',12);
title('垂向-校正边际谱','fontsize',12);