Matlab信号处理：频域分析中的包络谱
包络谱是旋转机械故障诊断中一种重要的分析手段。顾名思义，包络谱就是信号包络的频谱分析结果，它主要针对调幅信号的解调。通常，先对原始信号去均值，即去趋势化，采用希尔伯特变换，将原始信号转换为解析信号，然后取模转换为上包络信号，再次去均值，然后进一步采用傅里叶变换获得包络谱。
1.具体案例
包络谱是幅值调制信号分析的重要手段。幅值调制是将一个高频载波信号与被测信号（调制信号）相乘，使得高频信号的幅值随着被测信号的变化而变化。本文所分析的调制信号，如下图所示：
其中，x1为高频载波信号，x2为调制信号，也是被测信号，20为直流偏置量，采样频率为512，采样点数为2048，即采样时间为4s。原始信号的波形如下图所示：
从信号的时域波形能看出，信号的上包络近似为调制信号的绝对值波形。从信号的频域分析中能发现，99.5Hz和100.5Hz的调制频率，其为100Hz和0.5Hz调制的结果。此外，频域中在0Hz也存在显眼的直流分量。
对上述信号，进行去均值，然后进行希尔伯特变换，获得解析信号，再取模获得上包络信号，结果如下图左所示。对上包络信号去均值，然后进一步傅里叶变换的频域如下图所示。
从图中可发现，获得的上包络信号基本为原始信号的包络，达到预期要求。从上图中右图中能发现，包络谱中存在1Hz的基频，它实际上就是包络信号的频率，因为是包络信号是以余弦信号（被测信号）的半个周期为一个周期的，所以包络信号的频率为1Hz。
为了进一步验证包络谱分析正确，采用matlab的包络谱函数进行验证，envspectrum函数设置Method为hilbert（希尔伯特），分析带宽尽可能设置大，最终为[frequ(2)，frequ(end-1)]，frequ表示频域分析的横轴，带宽范围为最低频域分辨率到采样频率减去最低频域分辨率。采用matlab的包络谱函数envspectrum获得包络信号（图左）和包络谱（图右）如下图所示：
从上图中能发现，所获得包络信号、包络谱，和matlab自带函数获得包络信号和包络谱基本重合。具体局部细节放大如下图所示：
从定量结果来看，两者之间的差值也是很小的，可以忽略。综上分析，所以本文包络谱分析的过程可以看作是准确的。
除此之外，在包络谱分析的过程中，涉及到2个去均值的地方，一个在进行包络谱分析前，另外一个在最后的傅里叶变换前。在傅里叶变换时均值不去在包络谱中仅仅是一个直流分量，影响不大。但是，如果包络谱分析前不去掉均值、或未去到位，希尔伯特提取包络谱效果非常不佳，还会出现一些伪分量。
上图左图中能发现，由于未去均值，希尔伯特变换取模后的包络信号，和原始信号基本一致。从右图能看到，100Hz处的调制频率没有被完全解调，在200Hz处出现了原始信号中没有的伪分量，这干扰了信号的分析。由此可见，在希尔伯特变换前去均值或去趋势化的重要性，因为这可能导致信号的包络提取失败。
2.具体代码
主函数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,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上述是一个测试信号包络谱分析的实例，测试信号是一个幅值调制的实例，采用包络谱分析其调制信号的过程。同时采用matlab的envspectrum函数进行算法的验证，验证手动计算的正确性，此外通过实例说明了希尔伯特变换提取包络前去均值的重要性。