pengcheng5251 发表于 2011-4-18 22:18

求助:基于EMD的能量算子在MATLAB中的实现

我现在在做包络解调,最近学习了基于EMD的能量算子解调,但是不会在matlab中编写程序,求教程序以及原函数等,最好有仿真信号的例子。谢谢大家!

王晓华 发表于 2011-5-8 10:30

回复 1 # pengcheng5251 的帖子

正在研究,同求助!

impulsesignal 发表于 2011-11-4 15:52

楼主问题解决了吗?能否和大家分享一下哦{:{13}:}

金屋狗剩 发表于 2011-11-6 18:25

请求结果分享!学习中....

manguoyong 发表于 2011-11-15 14:40

请求结果分享!学习中....

石头王石头 发表于 2012-9-18 10:41

解调信号,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))

石头王石头 发表于 2012-9-18 10:42

这是我在一篇博客上看到的,分享一下,公式没复制上,不过不影响

zzdxhuangcj 发表于 2012-9-18 16:39

《机械故障诊断的Hilbert-huang变换方法》于德介、程军圣主编对基于EMD的能量算子解调方法介绍的比较详细,应该也发过类似的文章。

念家的猫 发表于 2012-10-26 21:29

这个貌似是黄09年提出来的新方法

shuihai707 发表于 2012-11-8 22:14

zzdxhuangcj 发表于 2012-9-18 16:39 static/image/common/back.gif
《机械故障诊断的Hilbert-huang变换方法》于德介、程军圣主编对基于EMD的能量算子解调方法介绍的比较详细, ...

有电子版的吗,能发一份吗?谢谢啦,1397245858@qq.com

华电机械 发表于 2013-7-10 13:03

石头王石头 发表于 2012-9-18 10:41 static/image/common/back.gif
解调信号,Hilbert是一种有效的方法,但它对信号由要求,需要满足Bedrosian Product Theorem,才不会有大的 ...

lz所说的hht-toolbox有么,,传我要一份呗,谢了,354773461@qq.com

hht123 发表于 2015-4-16 21:57

华电机械 发表于 2013-7-10 13:03
lz所说的hht-toolbox有么,,传我要一份呗,谢了,


这里面有

华电机械 发表于 2015-4-17 20:54

hht123 发表于 2015-4-16 21:57
这里面有

谢谢~~~~~~~~~~~~
页: [1]
查看完整版本: 求助:基于EMD的能量算子在MATLAB中的实现