关于EMD中mask函数的用法

zhangnan3509

苦思冥想，仍无结果，就此放弃，心有不甘！望有识之士能帮我解惑！100振币权作心意！

xray

% 估计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

eight

原帖由 zhangnan3509 于 2007-6-4 20:59 发表

                               
登录/注册后可看大图

苦思冥想，仍无结果，就此放弃，心有不甘！望有识之士能帮我解惑！100振币权作心意！

那篇应该是会议论文，所以很多东西没有讲清楚。最主要是集中在 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} 貌似跟时刻无关

xray

试着回答 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) 的计算公式，以及文章中五幅图的绘制，请各位高手指正。

xray

% 计算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));

xray

% 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-');

xray

% 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

xray

% 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-');

xray

% 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

xray

% 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--');

xray

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

eight

原帖由 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) 估计应该正确

xray

原帖由 eight 于 2007-6-20 10:19 发表

                               
登录/注册后可看大图

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

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

zhangnan3509

感谢xray,又给我们带了惊喜！

junzifei

模态混叠确实是个问题
不知道版主试过之后感觉这种方法怎么样？
感觉问题太多了！