查看: 5584|回复: 33

[HHT] 关于EMD中mask函数的用法

[复制链接]
发表于 2007-6-4 20:59 | 显示全部楼层 |阅读模式
100体能
苦思冥想,仍无结果,就此放弃,心有不甘!望有识之士能帮我解惑!100振币权作心意!

最佳答案

查看完整内容

% 估计MASK信号的频率 clear tic f1 = 1776; f2 = 1000; fs = 44100; n = 1:4410; x = sin(2*pi*f1*n/fs) + sin(2*pi*f2*n/fs); IMF = emd(x, 'MAXMODES', 3); lay = abs(fft(IMF(1,:))); [ pk pos ] = max(lay(1:end/2)); lay_f = (pos-1)*10; fm_est1 = xray_cal_f(IMF(1,:), fs); buf1 = sprintf('IMF(1,:)主要频率成分是%fHz, 估计得到的MASK信号频率是%fHz',lay_f, fm_est1); disp(buf1); ...
回复
分享到:

使用道具 举报

发表于 2007-6-4 20:59 | 显示全部楼层

xray_formular.m

% 估计MASK信号的频率

clear

tic
f1 = 1776;
f2 = 1000;
fs = 44100;
n = 1:4410;
x = sin(2*pi*f1*n/fs) + sin(2*pi*f2*n/fs);

IMF = emd(x, 'MAXMODES', 3);
        
lay = abs(fft(IMF(1,:)));
[ pk pos ] = max(lay(1:end/2));
lay_f = (pos-1)*10;   
fm_est1 = xray_cal_f(IMF(1,:), fs);
buf1 = sprintf('IMF(1,:)主要频率成分是%fHz, 估计得到的MASK信号频率是%fHz',lay_f, fm_est1);
disp(buf1);

lay = abs(fft(IMF(2,:)));
[ pk pos ] = max(lay(1:end/2));
lay_f = (pos-1)*10;   
fm_est2 = xray_cal_f(IMF(2,:), fs);
buf2 = sprintf('IMF(2,:)主要频率成分是%fHz, 估计得到的MASK信号频率是%fHz',lay_f, fm_est2);
disp(buf2)

xm1 = sin(2*pi*fm_est1*n/fs);
IMF = emd(x, 'MAXMODES', 3, 'MASK', xm1);
lay = abs(fft(IMF(1,:)));
[ pk pos ] = max(lay(1:end/2));
lay_f1 = (pos-1)*10;   
lay = abs(fft(IMF(2,:)));
[ pk pos ] = max(lay(1:end/2));
lay_f2 = (pos-1)*10;   
buf1 = sprintf('采用%fHz的MASK信号, 提取出来的IMF(1,:)是%fHz信号, IMF(2,:)是%fHz信号',fm_est1,lay_f1,lay_f2);
disp(buf1);
figure(1);hold on;
plot(IMF(1,:),'r');
plot(IMF(2,:),'b');

xm2 = sin(2*pi*fm_est2*n/fs);
IMF = emd(x, 'MAXMODES', 3, 'MASK', xm2);
lay = abs(fft(IMF(1,:)));
[ pk pos ] = max(lay(1:end/2));
lay_f1 = (pos-1)*10;   
lay = abs(fft(IMF(2,:)));
[ pk pos ] = max(lay(1:end/2));
lay_f2 = (pos-1)*10;   
buf2 = sprintf('采用%fHz的MASK信号, 提取出来的IMF(1,:)是%fHz信号, IMF(2,:)是%fHz信号',fm_est2,lay_f1,lay_f2);
disp(buf2);
figure(2);hold on;
plot(IMF(1,:),'r');
plot(IMF(2,:),'b');

disp('结论:MASK信号的频率选择,对IMF信号的提取有很大影响!');

toc
回复

使用道具 举报

发表于 2007-6-14 23:38 | 显示全部楼层


那篇应该是会议论文,所以很多东西没有讲清楚。最主要是集中在 mask 信号如何选取上。几个地方我都看不明白:

1) 4.1 节中,第二段的 the first unaltered IMF 是什么意思?
2) 能量加权平均(EWM)中的 k 如何取值?
3) An approximation of the frequency of the higher component can be found by finding the EWM of all points whose frequency is greater than \tilda{f}. 中的 of all points 是什么意思?每个时刻都有一个 EWM 吗?但是根据公式,\tilda{f} 貌似跟时刻无关
回复

使用道具 举报

发表于 2007-6-20 10:11 | 显示全部楼层

回复 #2 eight 的帖子

试着回答 eight 的几个问题:
1) the first unaltered IMF  在试验中发现,应该是指低频分量所对应的IMF;
2) 3) 能量加权平均(EWM)需要计算两次,第一次k=N (采样点总数), 计算得到f_1, 第二次选择所有大于f_1的频率f,把这些f和对应的a挑出来,再计算一次得到f_2。此时,f_2就是MASK信号的频率;

关于 “The use of a masking signal to improve empirical mode decomposition” 这篇文章, 我写了一些matlab的仿真代码,包括能量加权平均(EWM) 的计算公式,以及文章中五幅图的绘制,请各位高手指正。

评分

1

查看全部评分

回复

使用道具 举报

发表于 2007-6-20 10:12 | 显示全部楼层

xray_cal_f.m

% 计算mask信号的公式

function fm = cal_f(x, fs)

y = hilbert(x);
u = real(y);
v = imag(y);
a = sqrt(u.^2 + v.^2);

len = length(x);
f = zeros(1,len);
ts = 1/fs;
for k = 2:len-1
    f(k) = (u(k)*(v(k+1)-v(k-1))-v(k)*(u(k+1)-u(k-1)))/(2*pi*ts*(u(k)^2+v(k)^2));
end
f(1) = f(2);
f(len) = f(len-1);
f = fs * asin(f/fs);

fm_ = sum(a.*f.^2)/sum(a.*f);

index = find(f > fm_);
fm = sum(a(index).*f(index).^2)/sum(a(index).*f(index));
回复

使用道具 举报

发表于 2007-6-20 10:13 | 显示全部楼层

xray_fig1.m

% figure1
% 第3幅和第4幅图虽然趋势一致,但是尺度上大几个数量级

clear

tic
f1 = 1776;
f2 = 50:100:4000;
fs = 44100;
n = 1:4410;

len = length(f2);
lay1_f = zeros(1,len);
lay2_f = zeros(1,len);
lay1_am = zeros(1,len);
lay2_am = zeros(1,len);
lay1_av = zeros(1,len);
lay2_av = zeros(1,len);
err = zeros(1,len);

for k = 1:len
    f2(k)
    x = sin(2*pi*f1*n/fs) + sin(2*pi*f2(k)*n/fs);
    IMF = emd(x, 'MAXMODES', 3);
    if size(IMF,1) >= 1
        lay = abs(fft(IMF(1,:)));
        [ pk pos ] = max(lay(1:end/2));
        lay1_f(k) = (pos-1)*10;
        hht = abs(hilbert(IMF(1,:)));
        lay1_am(k) = mean(hht);
        lay1_av(k) = var(hht);
        err(k) = var(IMF(1,:)-x);
    end
    if size(IMF,1) >= 2
        lay = abs(fft(IMF(2,:)));
        [ pk pos ] = max(lay(1:end/2));
        lay2_f(k) = (pos-1)*10;   
        hht = abs(hilbert(IMF(2,:)));
        lay2_am(k) = mean(hht);
        lay2_av(k) = var(hht);
    end
end
toc
disp ok!

figure(1);
subplot(2,2,1);hold on;
plot(f2, lay1_f, 'r-');
plot(f2, lay2_f, 'b-.');
subplot(2,2,2);hold on;
plot(f2, lay1_am, 'r-');
plot(f2, lay2_am, 'b-.');
subplot(2,2,3);hold on;
plot(f2, lay1_av, 'r-');
plot(f2, lay2_av, 'b-.');
subplot(2,2,4);
plot(f2, err, 'r-');
回复

使用道具 举报

发表于 2007-6-20 10:13 | 显示全部楼层

xray_fig2.m

% figure2

clear

tic
f1 = 1776;
f2 = 1000;
fm = 2100;
fs = 44100;
n = 1:4410;
x = sin(2*pi*f1*n/fs);
x(1455:2955) = 0;
x = x + sin(2*pi*f2*n/fs);

IMF = emd(x, 'MAXMODES', 3);

figure(1);
subplot(5,1,1)
plot(n/fs, x, 'b-');
subplot(5,1,3)
plot(n/fs, IMF(1,:), 'b');
subplot(5,1,2)
plot(n/fs, IMF(2,:), 'r');

lay = abs(fft(IMF(1,:)));
[ pk pos ] = max(lay(1:end/2));
lay_f1 = (pos-1)*10;
lay = abs(fft(IMF(2,:)));
[ pk pos ] = max(lay(1:end/2));
lay_f2 = (pos-1)*10;
buf1 = sprintf('IMF(1,:)是%fHz信号, IMF(2,:)是%fHz信号',lay_f1,lay_f2);
disp(buf1);

IMF = emd(x,  'MAXMODES', 3, 'MASK', sin(2*pi*fm*n/fs));
subplot(5,1,5)
plot(n/fs, IMF(1,:), 'b');
subplot(5,1,4)
plot(n/fs, IMF(2,:), 'r');

lay = abs(fft(IMF(1,:)));
[ pk pos ] = max(lay(1:end/2));
lay_f1 = (pos-1)*10;
lay = abs(fft(IMF(2,:)));
[ pk pos ] = max(lay(1:end/2));
lay_f2 = (pos-1)*10;
buf1 = sprintf('IMF(1,:)是%fHz信号, IMF(2,:)是%fHz信号',lay_f1,lay_f2);
disp(buf1);

toc
回复

使用道具 举报

发表于 2007-6-20 10:14 | 显示全部楼层

xray_fig3.m

% figure3
% 按照前面提供的公式计算,MASK信号频率应该为4417.612358Hz
% 如果按文中的MASK信号频率,趋势不太像
% 如果按修正的MASK信号频率,曲线不太像
% 而且误差的数量级也有较大差别

clear

tic

f1 = 3060;
f2 = 2050;
fs = 44100;
n = 1:4410;
x = sin(2*pi*f1*n/fs) + sin(2*pi*f2*n/fs);
s = sin(2*pi*f1*n/fs);

fm = [ 3160 3360 3660 3860 ];
% fm = [ 4160 4360 4660 4860 ];
for k = 1:4
    xm(k,:) = sin(2*pi*fm(k)*n/fs);
end

a = 0.1:0.2:3;
for k = 1:length(a)
    a(k)
    IMF = emd(x, 'MAXMODES', 3);
    err_1(k) = var(IMF(1,:)-s);
    err_2(k) = var(IMF(2,:)-s);

    for m = 1:4
        IMF = emd(x, 'MAXMODES', 3, 'MASK', a(k)*xm(m,:));
        errm_1(m, k) = var(IMF(1,:)-s);
        errm_2(m, k) = var(IMF(2,:)-s);            
    end

end
toc
disp ok!

figure(1); hold on;
plot(a, err_1, 'r');
plot(a, errm_1(1,:), 'b:');
plot(a, errm_1(2,:), 'b-.');
plot(a, errm_1(3,:), 'b--');
plot(a, errm_1(4,:), 'b-');

figure(2); hold on;
plot(a, err_2, 'r');
plot(a, errm_2(1,:), 'b:');
plot(a, errm_2(2,:), 'b-.');
plot(a, errm_2(3,:), 'b--');
plot(a, errm_2(4,:), 'b-');
回复

使用道具 举报

发表于 2007-6-20 10:15 | 显示全部楼层

xray_fig4.m

% figure4
% 对于这幅图,MASK方法放大了其中一个分量的幅度,但是优势不如文章中的明显
% 而且文章中fm的选择,居然在0.5<f_1/f_2<2的范围内,连原始信号中最小频率的
% 2倍都不到,怀疑此处有误,将f2=1000或者fm=2500,似乎效果略好一些。

clear

tic

f1 = 1776;
f2 = 1200;
fm = 2100;
fs = 44100;
n = 1:3087;
x = sin(2*pi*f1*n/fs) + sin(2*pi*f2*n/fs);

IMF = emd(x,  'MAXMODES', 3);
figure(1);
subplot(5,1,1)
plot(n/fs, x, 'b-');
ylim([ -2 2 ])
subplot(5,1,3)
plot(n/fs, IMF(1,:), 'b');
ylim([ -2 2 ])
subplot(5,1,2)
plot(n/fs, IMF(2,:), 'r');
ylim([ -2 2 ])

IMF = emd(x,  'MAXMODES', 3, 'MASK', sin(2*pi*fm*n/fs));
subplot(5,1,5)
plot(n/fs, IMF(1,:), 'b');
ylim([ -2 2 ])
subplot(5,1,4)
plot(n/fs, IMF(2,:), 'r');
ylim([ -2 2 ])

toc
回复

使用道具 举报

发表于 2007-6-20 10:15 | 显示全部楼层

xray_fig5.m

% figure5
% 曲线趋势上差不多,但是数量级不对,而且MASK方法的优点不是特别明显

clear

tic
f1 = 1776;
f2 = 500:100:2500;
fs = 44100;
n = 1:4410;

fm1 = 2000;
fm2 = 2100;
xm1 = sin(2*pi*fm1*n/fs);
xm2 = sin(2*pi*fm2*n/fs);
s = sin(2*pi*f1*n/fs);

len = length(f2);
err_s = zeros(1,len);
err_x = zeros(1,len);
err_m1_s = zeros(1,len);
err_m1_x = zeros(1,len);
err_m2_s = zeros(1,len);
err_m2_x = zeros(1,len);
LAY = 1;
for k = 1:len
    f2(k)
    x = sin(2*pi*f1*n/fs) + sin(2*pi*f2(k)*n/fs);
    IMF = emd(x, 'MAXMODES', 3);
    IMF_M1 = emd(x, 'MAXMODES', 3, 'MASK', xm1);
    IMF_M2 = emd(x, 'MAXMODES', 3, 'MASK', xm2);
    if size(IMF,1) >= LAY
        err_s(k) = var(IMF(LAY,:)-s);
        err_x(k) = var(IMF(LAY,:)-x);
    end
    if size(IMF_M1,1) >= LAY
        err_m1_s(k) = var(IMF_M1(LAY,:)-s);
        err_m1_x(k) = var(IMF_M1(LAY,:)-x);
    end
    if size(IMF_M2,1) >= LAY
        err_m2_s(k) = var(IMF_M2(LAY,:)-s);
        err_m2_x(k) = var(IMF_M2(LAY,:)-x);
    end   
end
toc
disp ok!

figure(1);
subplot(2,1,1);hold on;
plot(f2, err_s, 'r-');
plot(f2, err_m1_s, 'b:');
plot(f2, err_m2_s, 'b--');
subplot(2,1,2);hold on;
plot(f2, err_x, 'r-');
plot(f2, err_m1_x, 'b:');
plot(f2, err_m2_x, 'b--');
回复

使用道具 举报

发表于 2007-6-20 10:18 | 显示全部楼层

小结

代码贴完了,其中有些结果和文章中的不太一致,标注也没有添加,权作抛砖引玉了,欢迎大家指正!

评分

1

查看全部评分

回复

使用道具 举报

发表于 2007-6-20 10:19 | 显示全部楼层
原帖由 xray 于 2007-6-20 10:11 发表
试着回答 eight 的几个问题:
1) the first unaltered IMF  在试验中发现,应该是指低频分量所对应的IMF;
2) 3) 能量加权平均(EWM)需要计算两次,第一次k=N (采样点总数), 计算得到f_1, 第二次选择所有大于 ...


1) 文中 the first unaltered IMF  是 y_1(t),即第一个 IMF,这是最高频的分量,这与低频分量是否矛盾?
2) 估计应该正确
回复

使用道具 举报

发表于 2007-6-20 10:24 | 显示全部楼层
原帖由 eight 于 2007-6-20 10:19 发表


1) 文中 the first unaltered IMF  是 y_1(t),即第一个 IMF,这是最高频的分量,这与低频分量是否矛盾?
2) 估计应该正确



如果用第一个IMF,即最高频的分量,算出来的结果和作者给出的MASK频率不相符合,如果有兴趣可以用我给出来的代码算来试试。
回复

使用道具 举报

 楼主| 发表于 2007-6-20 10:43 | 显示全部楼层

回复 #13 xray 的帖子

感谢xray,又给我们带了惊喜!
回复

使用道具 举报

发表于 2007-7-9 11:27 | 显示全部楼层
模态混叠确实是个问题
不知道版主试过之后感觉这种方法怎么样?
感觉问题太多了!
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表