基于G.Rilling的极值镜像延拓边界处理方法的改进
本帖最后由 529899778 于 2012-1-12 17:40 编辑之前看到《希尔伯特-黄变换方法的改进》中提出“平行延拓”的方法来对端点效应问题进行改善处理：
“利用端点处附近的两个相邻极值点（一极大值，以极小值）处斜率相等这一特性，认为在边端处定义出两个极值点，分别连接相邻的极大值与极小值，对包络线进行拟合”。
上图为第一个极值是极小值与最后一个极值是极大值的情况，其他三种情况类似。
针对图2中的情况，G.Rilling方法（1）若第一个采样值小于第一个极大值，则以第一个极小值点所在时间点镜像延拓； （2）若第一个采样值大于第一个极小值，则以第一个采样点所在时间点镜像延拓。本人欲结合两种方法，进行镜像延拓，例如情况（2）中，用平行延拓方法预测得到一个极小值点，镜像延拓后样条拟合再求均值。
通过仿真信号发现并没有很好的改善端点效应问题，想请教论坛的朋友们，是本人理解有错还是程序有问题。
程序附上 （后缀改为.m即可运行） 本帖最后由 529899778 于 2012-1-13 09:54 编辑
对程序稍微进行了修改，还是直接将相对G.Rilling的程序修改的部分直接贴出来。
% boundary conditions for interpolations :
%%%插值的边界条件
if indmax(1) < indmin(1)%第一个极值点是极大值
if m(1) > m(indmin(1))%如果第一个采样点大于第一极小值点，以第一个极大值为对称中心
lmax = fliplr(indmax(2:min(end,NBSYM+1)));
lmin = fliplr(indmin(1:min(end,NBSYM)));
lsym = indmax(1);
else%如果第一个采样点小于第一个极小值点，则认为第一个采样点为极小值点，以第一个采样点的预测点为对称中心
lmax = ;
lmin = ;
lsym = 1;
end
else%第一个极值点是极小值
if m(1) < m(indmax(1))%以第一个极小值为对称中心
lmax = fliplr(indmax(1:min(end,NBSYM)));
lmin = fliplr(indmin(2:min(end,NBSYM+1)));
lsym = indmin(1);
else%以第一个采样点为对称中心
lmax = ;
lmin = ;
lsym = 1;
end
end
%%%末尾序列与开头序列类似
if indmax(end) < indmin(end)%最后一个极值是极小值
if m(end) < m(indmax(end))%最后一个采样点小于最后一个极大值点，以最后一个极小值点为对称中心
rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
rmin = fliplr(indmin(max(end-NBSYM,1):end-1));
rsym = indmin(end);
else%最后一个采样点大于最后一个极值点，认为该点为极大值点，以该点的预测点为对称中心
rmax = ;
rmin = ;
rsym = lx;
end
else
if m(end) > m(indmin(end))
rmax = fliplr(indmax(max(end-NBSYM,1):end-1));
rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
rsym = indmax(end);
else
rmax = ;
rmin = ;
rsym = lx;
end
end
%%%将序列根据中心对称镜像延拓到两边
tlmin = 2*t(lsym)-t(lmin);
tlmax = 2*t(lsym)-t(lmax);
trmin = 2*t(rsym)-t(rmin);
trmax = 2*t(rsym)-t(rmax);
% in case symmetrized parts do not extend enough
%%%如果对称的部分没有足够的极值点
if tlmin(1) > t(1) || tlmax(1) > t(1)%对折后的序列没有超出原序列的范围
if lsym == indmax(1)
lmax = fliplr(indmax(1:min(end,NBSYM)));
else
lmin = fliplr(indmin(1:min(end,NBSYM)));
end
if lsym == 1%以第一个采样点为对称中心镜像，但对折后的序列不超出原序列的范围的这种情况不应该出现，若出现，则程序中止
error('bug')
end
lsym = 1;%直接以第一采样点为对称中心取镜像
tlmin = 2*t(lsym)-t(lmin);
tlmax = 2*t(lsym)-t(lmax);
end
%%%序列末尾情况与序列开头类似
if trmin(end) < t(lx) || trmax(end) < t(lx)
if rsym == indmax(end)
rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
else
rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
end
if rsym == lx
error('bug')
end
rsym = lx;
trmin = 2*t(rsym)-t(rmin);
trmax = 2*t(rsym)-t(rmax);
end
%%%延拓点上的取值
mlmax =m(lmax);
mlmin =m(lmin);
mrmax =m(rmax);
mrmin =m(rmin);
tmin = ;
tmax = ;
mmin = ;
mmax = ;
if indmax(1)<indmin(1)%第一个极值是极大值
if m(1)<m(indmin(1))%第一个采样值小于第一个极小值，需要预测一个极小值
slope=(m(indmin(1))-m(1))/(indmin(1)-1);
b=m(indmax(1))-slope*indmax(1);
y1_lpred=slope+b;
mmax(length(mlmax))=y1_lpred;
end
else第一个极值是极小值
if m(1)>m(indmax(1))%第一个采样值大于第一个极大值，需要预测一个极大值
slope=(m(indmax(1))-m(1))/(indmax(1)-1);
b=m(indmin(1))-slope*indmin(1);
y2_lpred=slope+b;
mmin(length(mlmin))=y2_lpred;
end
end
%%%序列末尾情况与开头情况类似
if indmax(end)<indmin(end)
if m(end)>m(indmax(end))
slope=(m(end)-m(indmax(end)))/(lx-indmax(end));
b=x(indmin(end))-slope*indmin(end);
y1_rpred=slope*lx+b;
mmin(length(mlmin)+length(m(indmin))+1)=y1_rpred;
end
else
if m(end)<m(indmin(end))
slope=(m(end)-x(indmin(end)))/(lx-indmin(end));
b=m(indmax(end))-slope*indmax(end);
y2_rpred=slope*lx+b;
mmax(length(mlmax)+length(m(indmax))+1)=y2_rpred;
end
end
完整程序
function = boundary_conditions(indmin,indmax,t,x,z,nbsym)
lx=length(x);
indmax=;
indmin=;
zlmax=x(1);
if length(indmax)>=4
slope1=(x(indmax(3))-x(indmax(2)))/(indmax(3)-indmax(2));
tmp1=x(indmax(2))-slope1*(indmax(2)-indmax(1));
if tmp1>x(1)
zlmax=tmp1;
end
end
zrmax=x(end);
if length(indmax)>=4
slope2=(x(indmax(end-1))-x(indmax(end-2)))/(indmax(end-1)-indmax(end-2));
tmp2=x(indmax(end-1))+slope2*(indmax(end)-indmax(end-1));
if tmp2>x(end)
zrmax=tmp2;
end
end
zlmin=x(1);
if length(indmin)>=4
slope3=(x(indmin(3))-x(indmin(2)))/(indmin(3)-indmin(2));
tmp3=x(indmin(2))-slope3*(indmin(2)-indmin(1));
if tmp3<x(1)
zlmin=tmp3;
end
end
zrmin=x(end);
if length(indmin)>=4
slope4=(x(indmin(end-1))-x(indmin(end-2)))/(indmin(end-1)-indmin(end-2));
tmp4=x(indmin(end-1))+slope4*(indmin(end)-indmin(end-1));
if tmp4<x(end)
zrmin=tmp4;
end
end
tmin = t(indmin);
tmax = t(indmax);
zmin = ;
zmax = ;
end 这是我根据Huang在《Ensemble Empirical Mode Decomposition: a Noise-Assisted Data Analysis Method》中给出的算法写的程序，感觉不错。发出来给大家用。 playfish 发表于 2012-2-4 12:45 static/image/common/back.gif
这是我根据Huang在《Ensemble Empirical Mode Decomposition: a Noise-Assisted Data Analysis Method》中给 ...
非常感谢 先把您的程序尽量读懂 再请教您 playfish 发表于 2012-2-4 12:45 static/image/common/back.gif
这是我根据Huang在《Ensemble Empirical Mode Decomposition: a Noise-Assisted Data Analysis Method》中给 ...
您提到的文章之前下载过，Huang的EEMD的程序我也有，程序基本看懂了，就是用邻近端点的两个极值点求取一个值来代替端点，想必《希尔伯特-黄变换方法的改进》的作者也是看过这篇文章的，由于时间较紧张，文章没有细读，想请教您这样处理端点的理论依据是什么，您用该该方法所得记过的效果如何？ 回复 4 # playfish 的帖子
我根据你的程序对G.Rilling的程序进行了修改，然后对下边的信号进行imf分解：
fs=1000;
f1=50;
f2=100;
t=1/fs:1/fs:2;
x1=2*sin(2*pi*f1*t+0.2*pi);
x2=5*cos(2*pi*f2*t+0.5*pi);
x=;
clear fs f1 f2 t x1 x2
分解得到5个imf分量，第一个imf分量基本与原信号一致，后4个imf分量都是低频的，但最后一个imf分量明显不是低频走势，因为不是单调的，越来越困惑了........... 529899778 发表于 2012-2-9 16:59 static/image/common/back.gif
您提到的文章之前下载过，Huang的EEMD的程序我也有，程序基本看懂了，就是用邻近端点的两个极值点 ...
更正下，用的是与端点邻近的四个极值点。 回复 8 # 529899778 的帖子
g rilling的程序现在看已经不能用了。最好自己写一个。目前最新的进展指出，停止条件就是迭代10次就停止最好。然后NHT分离AM和FM。再对FM做DQ求出瞬时频率。 回复 6 # 529899778 的帖子
比rilling的镜像法更快收敛。边界抖动差不多。 回复 9 # playfish 的帖子
想请问playfish哪篇文献中有对“停止条件就是迭代10次就停止最好”这个问题的说明，确实，对仿真信号进行分解程序运行时间还是较短的，如果对实测信号分解时间就比较长了，这个需要改进下。 回复 10 # playfish 的帖子
对，我现在最头疼的问题就是端点效应，3个月的学习和编程，用6种方法（晚点我会发帖把我做的写出来，希望playfish能帮帮我）都无法得到一个好点的效果，以至于毕设无法往下做，真是悲了个剧啊...... 回复 11 # 529899778 的帖子
注意看Huang的论文列表，2010年的两篇就是讨论这个的
http://rcada.ncu.edu.tw/research1_clip_reference.htm 回复 13 # playfish 的帖子
谢谢，确实，还是应该看看老Huang的文章，毕竟是他搞出来的算法，要权威些... playfish 发表于 2012-2-4 12:44 static/image/common/back.gif
function = boundary_conditions(indmin,indmax,t,x,z,nbsym)
lx=length(x);
...
这个貌似出不来结果啊，亲{:3_47:}
页:
[1]
2