信号处理——频域分析(包络谱)
信号处理——频域分析(包络谱)
包络谱是旋转机械故障诊断中一种重要的分析手段，主要用于调幅信号的解调分析。本文将详细介绍包络谱的原理、具体实现步骤以及注意事项，并通过一个具体的调幅信号案例进行演示。
1. 概述
包络谱是旋转机械故障诊断中一种重要的分析手段。顾名思义，包络谱就是信号包络的频谱分析结果，它主要针对调幅信号的解调。
本文绘制包络谱的方法的技术路线如下：
首先，先对原始信号去均值，也叫去趋势化（很重要，后面探讨了去均值和不取均值包络谱的差异），采用希尔伯特变换，将原始信号转换为解析信号，然后取模转换为上包络信号，再次去均值，然后进一步采用傅里叶变换获得包络谱。
该内容参考了一些资料：
1、书籍：MATLAB数字信号处理85个实用按比例——入门到进阶 宋知用 编著
2、书籍：机械工程测试技术基础 第3版 熊诗波编著
2、matlab官网的函数解释：
https://ww2.mathworks.cn/help/signal/ref/envspectrum.html?searchHighlight=envspectrum&s_tid=srchtitle_support_results_1_envspectrum
代码采用了Matlab 2024a进行运行，欢迎大家测试和提出问题！
2. 具体案例
包络谱是幅值调制信号分析的重要手段。幅值调制是将一个高频载波信号与被测信号（调制信号）相乘，使得高频信号的幅值随着被测信号的变化而变化。本文所分析的调制信号，如下图所示：
这里，x1为高频载波信号，x2为调制信号，也是被测信号，20为直流偏置量，采样频率为512，采样点数为2048，即采样时间为4s。原始信号的波形如下图所示：
从信号的时域波形能看出，信号的上包络近似为调制信号的绝对值波形。从信号的频域分析中能发现，99.5Hz和100.5Hz的调制频率，他为100Hz和0.5Hz调制的结果。此外，频域中在0Hz也存在显眼的直流分量。
对上述信号，进行去均值，然后进行希尔伯特变换，获得解析信号，再取模获得上包络信号，结果如下图所示。对上包络信号去均值，然后进一步傅里叶变换的频域如下图所示。
从上图中第一个子图中能发现，获得的上包络信号基本为原始信号的包络，达到预期要求。从上图中第二个子图中能发现，包络谱中存在1Hz的基频，它实际上就是包络信号的频率，因为是包络信号是以余弦信号（被测信号）的半个周期为一个周期的，所以包络信号的频率为1Hz。
因此本次包络谱分析的过程是正确的，也欢迎大家指出毛病哈，学艺不精，欢迎指正！
为了进一步验证我的包络谱分析正确，我采用matlab的包络谱函数envspectrum，进行验证，采用envspectrum函数设置Method为hilbert（希尔伯特）,分析带宽我尽可能设置大，最终为[frequ(2)，frequ(end-1)]，frequ表示频域分析的横轴，带宽范围为最低频域分辨率到采样频率减去最低频域分辨率，在envspectrum函数中，带宽无法设置为0到采样频率，会报错，没办法！也欢迎大家提供更好的验证方法。
采用matlab的包络谱函数envspectrum获得包络信号（第一个子图）和包络谱（第二个子图）如下图所示：
从上图中能发现，所获得包络信号、包络谱，和matlab自带函数获得包络信号和包络谱基本重合。具体局部细节放大如下图所示：
上述局部细节也能发现，手动计算方法和matlab自带函数的结果相差不多，下图为平均偏差结果：（平均偏差为两数差的绝对值的平均数）
从定量结果来看，两者之间的差值也是很小的，可以忽略。综上分析，所以本文包络谱分析的过程可以看作是准确的，欢迎大家提出批评指正！
除此之外，我还想和大家分享一下包络谱分析的注意事项。在包络谱分析的过程中，涉及到2个去均值的地方，一个在进行包络谱分析前，另外一个在最后的傅里叶变换前。在傅里叶变换时均值不去在包络谱中仅仅是一个直流分量，影响不大。但是，如果包络谱分析前不去掉均值、或未去到位，希尔伯特提取包络谱效果非常不佳，还会出现一些伪分量。下图为原始信号未去均值获得包络谱：
上图左图中能发现，由于未去均值，希尔伯特变换取模后的包络信号，和原始信号基本一致，下图为局部的一些放大图，两者相差的并不多。上图右图为左图包络信号进一步进行傅里叶变换的结果，没有去均值，所以信号中有直流分量，这不重要。从右图能看到，100Hz处的调制频率没有被完全的解调，在200Hz处出现了原始信号中没有的伪分量，这干扰了信号的分析。
由此可见，在希尔伯特变换前去均值或去趋势化的重要性，因为这个可能导致信号的包络提取失败。不过一般振动信号的偏置都接近0，不过最好还是去均值吧，严谨一点。
这个均值的发现要感谢书籍：MATLAB数字信号处理85个实用按比例——入门到进阶 宋知用 编著，要不然我也不会对去均值这么敏感。
3. 具体代码
代码主要分为两个：
1、envelope_main2.m（主函数，用于包络谱分析）
2、frequ_am_phase.m（幅值谱和相位谱计算函数，在本文中主要就是傅里叶变换使用）
说明：
frequ_am_phase.m函数
调用形式：
[freq,P1,Theta]=frequ_am_phase(y,fs,tol)
输入：
信号y，矩阵，行X列=单个信号的采样索引X信号数，比如信号的大小为8192X12，表示一个有12个信号的数据矩阵，每个信号长度为8192。注意，如果仅有一个信号，则y应该是一个列向量。同时，y的行数尽量为偶数，奇数的话会引起程序索引的警告。
fs:采样频率，标量
tol:相位阈值参数，标量（可缺省，下面讲解）
输出
freq:表示幅值谱的横轴，向量，HzX1
P1:表示幅值谱的纵轴，矩阵，单个信号的采样索引X信号数，类比信号y
Theta:表示相位谱的纵轴，矩阵，单个信号的采样索引X信号数，类比信号y
envelope_main2.m（主函数）
%% 信号的包络谱分析
%% 作者：冷漠
%% 时间：2024年7月10日
%% 关注公众号 ："故障诊断与寿命预测工具箱",每天进步一点点
%%
clc
clear all
close all
Fs = 512; % Sampling frequency 采样频率
T = 1/Fs; % Sampling period 采样周期
L = 2^11; % Length of signal 信号长度
t = (0:L-1)*T; % Time vector 时间向量
%% 仿真信号的构建
x1=cos(2*pi*100*t); %高频载波信号
x2=10.*cos(2*pi*0.5*t); %调制信号，也叫被测信号
y=(x1.*x2+20)'; %原始信号
[frequ,P1,~]=frequ_am_phase(y,Fs); %幅值谱分析
figure
subplot 211
plot(t,y,'b');title('信号时域分析'); axis tight;xlabel('t/s');ylabel('Amplitude')
subplot 212
plot(frequ,P1,'b');title('信号频域分析'); axis tight;xlabel('f/Hz');ylabel('Amplitude')
axis([0 max(frequ) min(P1) 6]);
%% 未去均值的包络谱分析
Hx= hilbert(y); %希尔伯特变换
f1=abs(Hx);
[frequ,P2,~]=frequ_am_phase(f1,Fs);
%作图
figure;
plot(t,f1,'r','LineWidth',1);hold on
plot(t,y,'b');
title('包络信号时域分析'); axis tight;xlabel('t/s');ylabel('Amplitude')
legend('未去均值获得的上包络信号','未去均值的原始信号')
figure
plot(frequ,P2,'b');title('信号的包络谱分析'); axis tight;xlabel('f/Hz');ylabel('Amplitude')
axis([0 max(frequ) 0 6]);
%% 去均值的包络谱分析
Hx1= hilbert(y-mean(y)); %希尔伯特变换
f2=abs(Hx1);
f3=f2-mean(f2);
[frequ,P3,~]=frequ_am_phase(f3,Fs);
%作图
figure;
plot(t,f2,'r','LineWidth',1);hold on;
plot(t,y-mean(y),'b');
title('包络信号时域分析'); axis tight;xlabel('t/s');ylabel('Amplitude');
legend('上包络信号','去均值的原始信号')
figure
plot(frequ,P3,'b');title('信号的包络谱分析'); axis tight;xlabel('f/Hz');ylabel('Amplitude')
axis([0 max(frequ) 0 6]);
%% Matlab自测函数
[es frequ_matlab env]=envspectrum(y,Fs,"Method","hilbert",'Band',[frequ(2) frequ(end-1)]);
%作图
figure;
plot(t,f3,'k','LineWidth',1);hold on;
plot(t,env,'r','LineWidth',1);
legend('envspectrum函数获得上包络信号','去均值的上包络信号')
title('包络信号时域对比'); axis tight;xlabel('t/s');ylabel('Amplitude');
fprintf('matlab的envspectrum函数获得包络与手写函数之间的差为%f \n',sum(abs(env-f3))/length(env));
figure
plot(frequ,P3,'k');hold on;
plot(frequ,es,'r');
legend('手动获得的包络谱','envspectrum函数获得的包络谱');
title('信号的包络谱对比'); axis tight;xlabel('f/Hz');ylabel('Amplitude')
axis([0 25 0 6]);
fprintf('matlab的envspectrum函数获得包络谱与手写函数之间的差为%f \n',sum(abs(es(2:end-1)-P3(2:end-1)))/length(es(2:end-1)));
frequ_am_phase.m（幅值谱和相位谱计算函数）
function [freq,P1,Theta]=frequ_am_phase(y,fs,tol)
%% 作者：冷漠
%% 时间：2024年6月01日
%% 关注公众号 ："故障诊断与寿命预测工具箱",每天进步一点点
%% 绘制信号频域的幅值谱和相位谱
%% 参数解释：
% y: 表示输入信号，它可以为一个矩阵，行X列,具体为单个信号的采样索引X信号数
% 比如y的大小为8192X12，表示一个有12个信号的数据矩阵，每个信号长度为8192
% 注意，如果仅有一个信号，则y应该是一个列向量
% 同时，y的行数尽量为偶数，奇数的话会引起程序索引的警告
% fs:表示采样频率
% tol:相位阈值参数
% freq:表示幅值谱的横轴
% P1:表示幅值谱的纵轴
% Theta:表示相位谱的纵轴
if nargin==2
tol=1e-6; %计算误差的默认阈值
end
L=size(y,1); % 信号长度
% Y=fft(y,2^nextpow2(L)); % FFT 快速傅里叶变换
Y=fft(y,L); % FFT 快速傅里叶变换
freq=(0:L/2)*fs/L; % 设置频率刻度 横轴Hz
%幅值谱
P2 = abs(Y/L);
P1 = P2(1:L/2+1,:);
P1(2:end-1,:) = 2*P1(2:end-1,:); %纵轴 幅值
%相位谱
P2(2:end-1,:)=2*P2(2:end-1,:);
for i=1:size(Y,2)
Y(P2(:,i)<tol,i) = 0;
theta(:,i) = angle(Y(:,i))/pi;
end
Theta=theta(1:L/2+1,:);
end
4. 细节说明
1、在采用希尔伯特变换提取信号包络时，要注意信号的包络是否自己想要的，他直接决定了包络谱的结果是否靠谱！一定要记得去均值或者去趋势化！很重要！重要！重要！
2、关于matlab自带函数envspectrum挺好用的，封装不错，大家可以试一试进行包络谱分析，但是需要设置带宽。按照代码中我这样设置是可以的，frequ就是频域分析的横轴坐标，代码如下：frequ=(0:L/2)*fs/L;L为信号长度，fs采用频率。
5. 总结
上述是一个测试信号包络谱分析的例子，测试信号是一个幅值调制的例子，采用包络谱分析其调制信号的过程。同时，采用matlab的envspectrum函数进行算法的验证，验证手动计算的正确性。此外，通过例子说明了希尔伯特变换提取包络前去均值的重要性。
包络谱分析的步骤：1去均值，2希尔伯特变换获得解析信号，3取模获得包络信号，4去均值进行傅里叶变换，获得包络谱。
6. 相关资料
附件
opgs
1、上述源码
①代码文件：
envelope_main2.m(主函数)；
frequ_am_phase.m（幅值谱和相位谱计算函数）
2、相关参考
①书籍：MATLAB数字信号处理85个实用按比例——入门到进阶 宋知用 编著
②书籍：机械工程测试技术基础 第3版 熊诗波编著
③matlab官网的函数解释：
更多内容，请关注公众号“故障诊断与寿命预测工具箱”，每天进步一点点。