关于EMD中mask函数的用法
苦思冥想，仍无结果，就此放弃，心有不甘！望有识之士能帮我解惑！100振币权作心意！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 原帖由 zhangnan3509 于 2007-6-4 20:59 发表 http://www.chinavib.com/forum/images/common/back.gif
苦思冥想，仍无结果，就此放弃，心有不甘！望有识之士能帮我解惑！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} 貌似跟时刻无关
回复 #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) 的计算公式，以及文章中五幅图的绘制，请各位高手指正。
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));
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-');
xray_fig2.m
% figure2clear
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_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-');
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
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--');
小结
代码贴完了，其中有些结果和文章中的不太一致，标注也没有添加，权作抛砖引玉了，欢迎大家指正！ 原帖由 xray 于 2007-6-20 10:11 发表 http://www.chinavib.com/forum/images/common/back.gif试着回答 eight 的几个问题：
1) the first unaltered IMF在试验中发现，应该是指低频分量所对应的IMF；
2) 3) 能量加权平均(EWM)需要计算两次，第一次k=N (采样点总数), 计算得到f_1， 第二次选择所有大于 ...
1) 文中 the first unaltered IMF是 y_1(t)，即第一个 IMF，这是最高频的分量，这与低频分量是否矛盾？
2) 估计应该正确 原帖由 eight 于 2007-6-20 10:19 发表 http://www.chinavib.com/forum/images/common/back.gif
1) 文中 the first unaltered IMF是 y_1(t)，即第一个 IMF，这是最高频的分量，这与低频分量是否矛盾？
2) 估计应该正确
如果用第一个IMF，即最高频的分量，算出来的结果和作者给出的MASK频率不相符合，如果有兴趣可以用我给出来的代码算来试试。
回复 #13 xray 的帖子
感谢xray,又给我们带了惊喜！ 模态混叠确实是个问题不知道版主试过之后感觉这种方法怎么样？
感觉问题太多了！