信号处理中的包络谱分析：原理、实现与注意事项

信号处理中的包络谱分析：原理、实现与注意事项

引用

CSDN

1.

https://m.blog.csdn.net/ynn4818172/article/details/140423173

1.概述

包络谱是旋转机械故障诊断中一种重要的分析手段。顾名思义，包络谱就是信号包络的频谱分析结果，它主要针对调幅信号的解调。

本文将介绍绘制包络谱的具体方法和技术路线：

1. 首先对原始信号进行去均值处理（也称去趋势化），这一步非常重要，后面会探讨其必要性。
2. 采用希尔伯特变换将原始信号转换为解析信号。
3. 取模运算获得上包络信号，并再次进行去均值处理。
4. 最后通过傅里叶变换获得包络谱。

该内容参考了以下资料：

1. 书籍：《MATLAB数字信号处理85个实用案例——入门到进阶》 宋知用 编著
2. 书籍：《机械工程测试技术基础》 第3版 熊诗波编著
3. MATLAB官网函数解释：envspectrum

代码采用MATLAB 2024a进行运行，欢迎读者测试和提出问题。

2.具体案例

包络谱是幅值调制信号分析的重要手段。幅值调制是将一个高频载波信号与被测信号（调制信号）相乘，使得高频信号的幅值随着被测信号的变化而变化。

本文分析的调制信号如下：

• x1为高频载波信号
• x2为调制信号（被测信号）
• 直流偏置量为20
• 采样频率为512Hz
• 采样点数为2048（即采样时间为4秒）

原始信号的波形如下：

从信号的时域波形可以看出，信号的上包络近似为调制信号的绝对值波形。在频域分析中可以发现99.5Hz和100.5Hz的调制频率，这是100Hz和0.5Hz调制的结果。此外，频域中在0Hz处存在显眼的直流分量。

对上述信号进行处理：

1. 去均值
2. 希尔伯特变换获得解析信号
3. 取模获得上包络信号

处理结果如下：

从上图第一个子图可以看出，获得的上包络信号基本符合预期。从第二个子图可以看出，包络谱中存在1Hz的基频，这是由于包络信号以余弦信号（被测信号）的半个周期为一个周期，因此包络信号的频率为1Hz。

为了验证分析的正确性，采用MATLAB的包络谱函数envspectrum进行验证，设置Method为hilbert（希尔伯特），分析带宽设置为最低频域分辨率到采样频率减去最低频域分辨率。

采用MATLAB的包络谱函数envspectrum获得的包络信号和包络谱如下：

从上图可以看出，手动计算方法和MATLAB自带函数的结果基本一致。局部细节放大图进一步证实了这一点：

定量分析显示，两者之间的平均偏差很小，可以忽略。因此，本文的包络谱分析过程可以认为是准确的。

包络谱分析注意事项

在包络谱分析过程中涉及两个关键的去均值步骤：

1. 在进行包络谱分析前
2. 在最后的傅里叶变换前

在傅里叶变换时，如果不去除均值，只会导致包络谱中出现一个直流分量，影响不大。但是，如果在包络谱分析前不去掉均值或未去到位，希尔伯特变换提取包络的效果会非常差，还可能出现伪分量。

下图展示了原始信号未去均值时获得的包络谱：

从上图左图可以看出，由于未去均值，希尔伯特变换取模后的包络信号与原始信号基本一致。从右图可以看出，100Hz处的调制频率没有被完全解调，在200Hz处出现了原始信号中没有的伪分量，这干扰了信号的分析。

由此可见，在希尔伯特变换前进行去均值或去趋势化的重要性，因为这可能导致信号的包络提取失败。虽然一般振动信号的偏置都接近0，但最好还是进行去均值处理，以保持严谨性。

这个发现主要得益于《MATLAB数字信号处理85个实用案例——入门到进阶》一书的启发。

3.具体代码

代码分为两个部分：

1. envelope_main2.m（主函数，用于包络谱分析）
2. frequ_am_phase.m（幅值谱和相位谱计算函数，在本文中主要用于傅里叶变换）

envelope_main2.m（主函数）

%% 信号的包络谱分析
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)
%% 绘制信号频域的幅值谱和相位谱
%% 参数解释： 
%     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,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=(0:L/2)*fs/L，其中L为信号长度，fs为采样频率。

5.总结

本文通过一个测试信号包络谱分析的例子，展示了如何分析调幅信号的调制信号。通过与MATLAB自带的envspectrum函数结果进行对比，验证了手动计算方法的正确性。此外，通过实例说明了希尔伯特变换提取包络前去均值的重要性。

包络谱分析的基本步骤是：

1. 去均值
2. 希尔伯特变换获得解析信号
3. 取模获得包络信号
4. 去均值后进行傅里叶变换，获得包络谱。