信号处理——频域分析(包络谱)

最新推荐文章于 2024-11-16 18:19:00 发布

在职研究僧

最新推荐文章于 2024-11-16 18:19:00 发布

阅读量6k 信号处理算法

本文链接：https://blog.csdn.net/ynn4818172/article/details/140423173

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