求助：基于EMD的能量算子在MATLAB中的实现
我现在在做包络解调，最近学习了基于EMD的能量算子解调，但是不会在matlab中编写程序，求教程序以及原函数等，最好有仿真信号的例子。谢谢大家！ 回复 1 # pengcheng5251 的帖子正在研究，同求助！ 楼主问题解决了吗？能否和大家分享一下哦{:{13}:} 请求结果分享！学习中.... 请求结果分享！学习中.... 解调信号，Hilbert是一种有效的方法，但它对信号由要求，需要满足Bedrosian Product Theorem，才不会有大的误差，而且Hilbert变换的计算过程中涉及加窗，也会给结果带来误差。最明显的效应是在非整周期采样情况下，Hilbert变换会带来边缘飞翼，边缘误差很大。
实际上，实现解调的数学方法是很多的，如能量算子、Mandelstam法、Shekel法、Prony法等，而其中的能量算子法简单使用，能有效克服Hilbert解调的误差。
能量算子法的定义如下
其离散算法主要有两种：DESA-1和DESA-2，它们的计算公式分别如下
DESA-1:
DESA-2:
其中的ψ为Teager算子
知道了这些，既可以编写程序实现瞬时频率与包络的计算了，当然，成熟的算法需要考虑各种细节，如边缘的延拓等，DESA-1的一个简单代码实现如下
% Discrete Energy Separation Algorithm 1(DESA-1)for estimation purpose on an AM signal.
clc
clear
close all
% 原始信号
k=0.8;
w=2*pi/256;
q=2*pi/12;
n = 1:512;
A = 1+k*cos(w*n);
x = A.*cos(q*n);
x = x(:); % 转化为列向量
% 能量算子计算瞬时频率与包络
y = circshift(x,-1)-x;
ey = y.^2-circshift(y,1).*circshift(y,-1); % Teager算子
ex = x.^2-circshift(x,1).*circshift(x,-1); % Teager算子
w = acos(1-(ey+circshift(ey,1))./(4*ex)); % 瞬时频率
a = sqrt(ex./(1-(1-(ey+circshift(ey,1))./(4*ex)).^2)); % 包络
figure
plot()
legend({'原始信号','瞬时频率','包络'})
可以看到，包络和瞬时频率都能正确的提取出来，只是没有顾及边缘效应，边端效果较差，但也只是差在几个点上。
网上查代码，还看到一个三点对称差分能量算子法，其算法思想如下
Phid.m ：
DEO3S.m：
三点对称差分：
求解该差分使用了传递函数方法：
系统：
其传递函数表示为：
将
作为该系统的输入，输出即为三点对称差分后的能量算子分离结果，即DESA中的：
DESA.m ：
对其代码稍作整理后，封装成如下函数
function = DESA(x)
% 三点对称差分能量算子求解包络
% 算法详细描述：
% 2Ψ
% |a(n)|=-------------------------------------
% sqrt( Ψ[ x(n+1)-x(n-1) ] )
phix = DEO3S(x);
phix1 = DEO3S(circshift(x',-1)'-circshift(x',1)');
U = 2*phix;
L = abs(sqrt(phix1));
A=abs(U./L); % 应去除分母为0的点
A(A==inf)=0;
pha = 1-phix1/2./phix;
pha(pha>1 & pha<-1) = 0;
f = 0.5*acos(pha);
end
function f=DEO3S(x)
% 使用传递函数法求解三点对称差分能量算子
% x：行向量
% H(z)=z(1+2*z^-1+z^-2)/4;
Px=Phid(x);
Ns=length(Px);
w=2*pi*(-Ns/2:Ns/2)/Ns;
w=; % 去0点
z=exp(1i*w);
Hz=z.*(1+2*z.^-1+z.^-2)/4;
Xz=fft(Px,Ns);
Xz=;%重新排列
Yw=Hz.*Xz;
Yw=;
f=real(ifft(Yw));
end
function f=Phid(x)
% Teager能量算子：y=x(n)^2-x(n-1)*x(n+1)
x=x';
f=x.^2-circshift(x,1).*circshift(x,-1);
f=f';
end
对此函数的测试如下
% function testPDESA
clc
clear
close all
fs=5120;
Ns=2048;
t=(0:Ns-1)/fs;
x2=exp(-5*t).*cos(2*pi*20*t); % 实际包络
x=x2.*sin(2*pi*500*t);
=DESA(x); % 能量算子解包络
yhilbert=abs(hilbert(x)); % Hilbert解包络
figure
subplot(2,1,1),
plot(t,)
title('解包络')
legend('能量算子包络','Hilbert包络','原始数据');
grid on
subplot(2,1,2),
plot(t,)
grid on
axis();
legend('能量算子包络误差','Hilbert包络误差');
title('解包络误差')
figure
plot(abs(f)/2/pi*fs)
title('能量算子解调得到的频率')
可以看到它求解得到的包络和Hilbert得到的包络相差无几，也可以看到能量算子法在一些过零点处有些误差，但还可以容忍，而Hilbert法的边端飞翼就确实太大了（能量算子只在边界一两个点上误差很大，毕竟能量算子中涉及离散差分，差分一次就少一点的信息，故造成边缘点的误差）；另外能量算子得到的瞬时频率貌似不照，其实如果用优化算法优化后，效果是可以得到很大提升的。
实际上，能量算子的实际应用可以利用一个matlab工具箱：hht_toolbox，它是一个希尔伯特黄变换的一个工具箱（不过貌似用的人不多，不如G.Rilling的代码好），里面有计算4个能量算子的函数：
desa：原始能量算子法
desa1：DESA-1法
desa1m：DESA-1法的平滑结果
desa2：DESA-2法
代码效果测试如下
% 测试各种能量算子算法的效果
clc
clear all
close all
ts = 0.001;
N = 200;
fs = 1/ts;
k = 0.8;
t = (0:N-1)*ts;
% 原始信号
fm = 10;
fc = 100;
Am = 1 + k*cos(2*pi*fm*t); % 包络
Am = 2+exp(-10*t).*cos(2*pi*fm*t); % 包络
Ac = cos(2*pi*fc*t);
y = Am.*Ac; % 信号调制
% 包络分析
yh = hilbert(y);
Ah = abs(yh); % 包络的绝对值
Ph = unwrap(angle(yh)); % 包络的相位
fh = diff(Ph)/2/pi; % 包络的瞬时频率，差分代替微分计算
% 能量算子解调
xn = y';
= desa(xn, ts);
= desa1(xn, ts);
= desa1m(xn, ts);
= desa2(xn, ts);
figure
subplot(211)
plot(t,)
title('能量算子解调振幅')
legend({'Hilbert包络','desa包络','desa1包络','desa1m包络','desa2包络','原始信号'})
fh(end+1) = fh(end);
subplot(212)
plot(t,)
title('能量算子解调频率')
legend({'Hilbert频率','desa频率','desa1频率','desa1m频率','desa2频率'})
% ylim()
figure
plot(t,)
title('能量算子解调振幅的误差')
legend({'Hilbert包络','desa包络','desa1包络','desa1m包络','desa2包络'})
sum(sum((A1-Am').^2))
sum(sum((A1m-Am').^2))
sum(sum((A2-Am').^2))
这是我在一篇博客上看到的，分享一下，公式没复制上，不过不影响 《机械故障诊断的Hilbert-huang变换方法》于德介、程军圣主编对基于EMD的能量算子解调方法介绍的比较详细，应该也发过类似的文章。 这个貌似是黄09年提出来的新方法 zzdxhuangcj 发表于 2012-9-18 16:39 static/image/common/back.gif
《机械故障诊断的Hilbert-huang变换方法》于德介、程军圣主编对基于EMD的能量算子解调方法介绍的比较详细， ...
有电子版的吗，能发一份吗？谢谢啦,1397245858@qq.com 石头王石头 发表于 2012-9-18 10:41 static/image/common/back.gif
解调信号，Hilbert是一种有效的方法，但它对信号由要求，需要满足Bedrosian Product Theorem，才不会有大的 ...
lz所说的hht-toolbox有么，，传我要一份呗，谢了，354773461@qq.com 华电机械 发表于 2013-7-10 13:03
lz所说的hht-toolbox有么，，传我要一份呗，谢了，
这里面有 hht123 发表于 2015-4-16 21:57
这里面有
谢谢~~~~~~~~~~~~
页:
[1]