极值波延拓matlab函数,[HHT] 基于G.Rilling的极值镜像延拓边界处理方法的改进

之前看到《希尔伯特-黄变换方法的改进》中提出“平行延拓”的方法来对端点效应问题进行改善处理：

“利用端点处附近的两个相邻极值点(一极大值，以极小值)处斜率相等这一特性，认为在边端处定义出两个极值点，分别连接相邻的极大值与极小值，对包络线进行拟合”。(在线时间不够，图片只能上传附件)

上图为第一个极值是极小值与最后一个极值是极大值的情况，其他三种情况类似。

针对图2中的情况，G.Rilling方法(1)若第一个采样值小于第一个极大值，则以第一个极小值点所在时间点镜像延拓；      (2)若第一个采样值大于第一个极小值，则以第一个采样点所在时间点镜像延拓。本人欲结合两种方法，进行镜像延拓，例如情况(2)中，用平行延拓方法预测得到一个极小值点，镜像延拓后样条拟合再求均值。

通过仿真信号发现并没有很好的改善端点效应问题，想请教论坛的朋友们，是本人理解有错还是程序有问题。

相对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 = [fliplr(indmax(1:min(end,NBSYM-1))),1];

lmin = [fliplr(indmin(1:min(end,NBSYM-1))),1];

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 = [fliplr(indmax(1:min(end,NBSYM-1))),1];

lmin = [fliplr(indmin(1:min(end,NBSYM-1))),1];

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 = [lx,fliplr(indmax(max(end-NBSYM+2,1):end))];

rmin = [lx,fliplr(indmin(max(end-NBSYM+2,1):end))];

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 = [lx,fliplr(indmax(max(end-NBSYM+2,1):end))];

rmin = [lx,fliplr(indmin(max(end-NBSYM+2,1):end))];

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 = [tlmin t(indmin) trmin];

tmax = [tlmax t(indmax) trmax];

mmin = [mlmin m(indmin) mrmin];

mmax = [mlmax m(indmax) mrmax];

if indmax(1)

if m(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)

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)

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

完整程序下载：

[本帖最后由 grace19861013 于 2012-1-13 10:10 编辑]