算法描述
1.找到所有极大值点。
2.找到极值点中最大的。
3.以该最大点为起点,向两边找,
计算当前点与右边所有点连线的斜率,取斜率最大的连线的端点作为包络线的下一个端点,同时更新当前点为该端点,重复,直到最后一个点。
计算当前点与左边所有点连线的斜率取斜率最小的连线的端点作为包络线的下一个端点,同时更新当前点为该端点,重复此过程,直到第一个点。
4.直到所有极大值点都加入包络线为止
5.得到包络线后用原曲线除以包络线得到包络线去除后的曲线。
运行结果
可见包络线提取和包络线去除的效果较好。
MATLAB代码
在main.m中调用Envelope函数生成包络线。
main.m
%Copyright (C) China University of Geoscience
%20181003670 Fu Zhixiang
%
%% 创建并绘制光谱曲线
data=[6,10,18,13,10,11,30,40,39,39,38,38,...
39.5,39,38,39,35,30,25,18,23,20,25,30,31,...
32,28,25,20,15,20,12,13,15,16,20,18,15,12,10,8,5,2,0,0]./100;
x=0.4:0.05:2.6;
plot(x,data);
hold on;
%% 创建包络线
[result_x,result_y]=Envelope(x,data);%调用自定义的提取包络线函数
%plot(result_x,result_y,'-+');
title('包络线');
%% 包络线插值
y=interp1(result_x,result_y,x,'pchip');
plot(x,y);
%原数据除以包络线
y1=data./y;
plot(x,y1);
Envelope.m
function [ output_x,output_y ] = Envelope( input_x,input_y )
%包络线提取
%Copyright (C) China University of Geoscience
%064181 FuZhixiang
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%找极大值索引
LocMaxInd=find(diff(sign(diff(input_y)))<0)+1;%极大值在input_x中的索引
LocMax=input_y(LocMaxInd);%存放所有极大值
len=length(LocMax);
input_len=length(input_x);
%找极大值中的最大值在LocMax中的索引
[~,MaxInd]=max(LocMax);
output_x=LocMaxInd(MaxInd);%该数组记录包络线端点在input中的索引
output_y=LocMax(MaxInd);%该数组记录选取构成包络线的极大值
%计算最大值右边的包络线
if(MaxInd~=len)
i=MaxInd;
while(i~=len)
k_arr=zeros(1,len-i);%k_arr的个数等于当前点右边极大值点的个数
%依次计算当前点右边所有点与当前点的连线,并存入k_arr
for j=(i+1):len
k_arr(j-i)=(LocMax(j)-LocMax(i))/(LocMaxInd(j)-LocMaxInd(i));
end
%更新i
[~,ind]=max(k_arr);
i=ind+i;
%将选出的端点加入结果
output_x=[output_x,LocMaxInd(i)];
output_y=[output_y,LocMax(i)];
end
end
%将末端的非极值点加入包络线
for j=(LocMaxInd(i)+1):input_len
output_x=[output_x,j];
output_y=[output_y,input_y(j)];
end
%计算最大值左边的包络线
if(MaxInd~=1)
i=MaxInd;
while(i~=1)
k_arr=zeros(1,i-1);%
%依次计算当前点右边所有点与当前点的连线,并存入k_arr
for j=1:(i-1)
k_arr(j)=(LocMax(j)-LocMax(i))/(LocMaxInd(j)-LocMaxInd(i));
end
[~,i]=min(k_arr);%更新i
%将选出的端点加入结果
output_x=[LocMaxInd(i),output_x];
output_y=[LocMax(i),output_y];
end
end
%将起始的非极值点加入包络线
for j=(LocMaxInd(i)-1):-1:1
output_x=[j,output_x];
output_y=[input_y(j),output_y];
end
output_x=input_x(output_x);
end