第1页 / 共10页
第2页 / 共10页
第3页 / 共10页
第4页 / 共10页
第5页 / 共10页
第6页 / 共10页
第7页 / 共10页
第8页 / 共10页
局部均值分解法.doc
function[PF,A,SI]=lmd1(m)
%最后一个 PF 是残余分量
%A 是瞬时赋值
%SI 是纯调频函数，求它的瞬时频率就是需要的频率
c=m;
k=0;
wucha1=0.001;
n_l=nengliang(m);
while 1
k=k+1;
a=1;
h=c;
[pf,a,si]=zhaochun1(a,h,wucha1);
c=c-pf;
PF(k,:)=pf;
A(k,:)=a;
SI(k,:)=si;
c_pos=pos(c);
n_c=nengliang(c);
n_pf=nengliang(pf);
%停止调节
%1.emd 用的是三次样条求包络，要求至少 3 个极值点，所以这里 c 的极值点个数也应
该至少为 3
%2.如果上一个 PF 的极值点数比下一个 PF 的极值点数少，说明结果也不正确（这个也
可以作为停止条件考虑进去）
%上面一句是否可以等价于当前 PF 的极值点个数一定要大于等于残量（c）的极值点个
数（目前是用这个作为停止条件的一个参考写入程序）
%3.当前 PF 分量的能量应该大于残量 c 的能量(这个有待商榷)
%4.残余能量不能大于信号能量
if n_pfn_l
PF(k+1,:)=c;
break
end
end
end
%%%%%%%%%%%%%%%%%%nengliang 函数%%%%%%%%%%
function p=nengliang(y)
% my=mean(y);
% p=(y-my).*(y-my);
% p=sum(p);
p=sum(abs(y).^2);
end
%%%%%%%%%%%%%找纯函数%%%%%%%%%%%%%%%%%
function [pf,a,si]=zhaochun1(a,h,wucha1)
chun_num=0;
while 1
chun_num=chun_num+1;
t=1:length(h);
[envmin,envmax,envmoy,indmin,indmax,indzer] = envelope(t,h,'spline');
mi=(envmax+envmin)./2;
ai=abs(envmax-envmin)./2;
a=a.*ai;
si=(h-mi)./ai;
h=si;
ai_funmax=max(ai);
ai_funmin=min(ai);
if (ai_funmax<=1+wucha1&&ai_funmin>=1-wucha1)
break
end
end
pf=a.*si;
chun_num
function [envmin, envmax,envmoy,indmin,indmax,indzer] = envelope(t,x,INTERP)
%computes envelopes and mean with various interpolations
NBSYM = 2;
% 边界延拓点数
DEF_INTERP = 'spline';
if nargin < 2
x = t;
t = 1:length(x);
INTERP = DEF_INTERP;
end
if nargin == 2
if ischar(x)
INTERP = x;
x = t;
t = 1:length(x);
end
end
if ~ischar(INTERP)
error('interp parameter must be ''linear'''', ''cubic'' or ''spline''')
end
if ~any(strcmpi(INTERP,{'linear','cubic','spline'}))
error('interp parameter must be ''linear'''', ''cubic'' or ''spline''')
end
if min([size(x),size(t)]) > 1
error('x and t must be vectors')
end
s = size(x);
if s(1) > 1
x = x';
end
s = size(t);
if s(1) > 1
t = t';
end
if length(t) ~= length(x)
error('x and t must have the same length')
end
lx = length(x);
[indmin,indmax,indzer] = extr(x,t);
%boundary conditions for interpolation
[tmin,tmax,xmin,xmax] = boundary_conditions(indmin,indmax,t,x,NBSYM);
% definition of envelopes from interpolation
envmax = interp1(tmax,xmax,t,INTERP);
envmin = interp1(tmin,xmin,t,INTERP);
if nargout > 2
envmoy = (envmax + envmin)/2;
end
end
function [tmin,tmax,xmin,xmax] = boundary_conditions(indmin,indmax,t,x,nbsym)
% computes the boundary conditions for interpolation (mainly mirror symmetry)
lx = length(x);
% 判断极值点个数
if (length(indmin) + length(indmax) < 3)
error('not enough extrema')
end
% 插值的边界条件
if indmax(1) < indmin(1)% 第一个极值点是极大值
if x(1) > x(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)));
lmin = [fliplr(indmin(1:min(end,nbsym-1))),1];
lsym = 1;
end
else
if x(1) < x(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)));
lsym = 1;
end
end
% 序列末尾情况与序列开头类似
if indmax(end) < indmin(end)
if x(end) < x(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 = fliplr(indmin(max(end-nbsym+1,1):end));
rsym = lx;
end
else
if x(end) > x(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 = fliplr(indmax(max(end-nbsym+1,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
% 延拓点上的取值
xlmax =x(lmax);
xlmin =x(lmin);
xrmax =x(rmax);
xrmin =x(rmin);
% 完成延拓
tmin = [tlmin t(indmin) trmin];
tmax = [tlmax t(indmax) trmax];
xmin = [xlmin x(indmin) xrmin];
xmax = [xlmax x(indmax) xrmax];
end
%------------------------------------------------------------------------------
---------------------
% 极值点和过零点位置提取
function [indmin, indmax, indzer] = extr(x,t);
%extracts the indices corresponding to extrema
if(nargin==1)
t=1:length(x);
end
m = length(x);
if nargout > 2
x1=x(1:m-1);
x2=x(2:m);
indzer = find(x1.*x2<0);
if any(x == 0)
iz = find( x==0 );
indz = [];
if any(diff(iz)==1)
zer = x == 0;
dz = diff([0 zer 0]);
debz = find(dz == 1);
finz = find(dz == -1)-1;
indz = round((debz+finz)/2);
else
indz = iz;
end
indzer = sort([indzer indz]);
end
end
d = diff(x);
n = length(d);
d1 = d(1:n-1);
d2 = d(2:n);
indmin = find(d1.*d2<0 & d1<0)+1;
indmax = find(d1.*d2<0 & d1>0)+1;
% when two or more consecutive points have the same value we consider only one
extremum in the middle of the constant area
% 当连续多个采样值相同时，把最中间的一个值作为极值点，处理方式与连 0 类似
if any(d==0)
imax = [];
imin = [];
bad = (d==0);
dd = diff([0 bad 0]);
debs = find(dd == 1);
fins = find(dd == -1);
if debs(1) == 1
if length(debs) > 1
debs = debs(2:end);
fins = fins(2:end);
else
debs = [];
fins = [];
end
end
if length(debs) > 0
if fins(end) == m
if length(debs) > 1
debs = debs(1:(end-1));
fins = fins(1:(end-1));
else
debs = [];
fins = [];
end
end
end
lc = length(debs);
if lc > 0
for k = 1:lc
if d(debs(k)-1) > 0
if d(fins(k)) < 0
imax = [imax round((fins(k)+debs(k))/2)];
end
else
if d(fins(k)) > 0
imin = [imin round((fins(k)+debs(k))/2)];
end
end
end
end
if length(imax) > 0
indmax = sort([indmax imax]);
end
if length(imin) > 0
indmin = sort([indmin imin]);
end
end
end
end
%%%%%%%%%%%%%%%%%%pos%%%%%%%%%%%%
function poss=pos(y)
%一个输出
%功能：找序列极值点位置坐标
%y:输入序列
%pos:极值点的序列位置坐标
[indmin,indmax]=position(y);
minmax=cat(2,indmin,indmax);
poss=sort(minmax);
end
%%%%%%%%%%%%%%position%%%%%%%
%返还极大值和极小值位置下标
%两个输出
function [indmin,indmax]=position(y)
m = length(y);
d = diff(y);
n = length(d);
d1 = d(1:n-1);
d2 = d(2:n);
indmin = find(d1.*d2<0 & d1<0)+1;
indmax = find(d1.*d2<0 & d1>0)+1;
if any(d==0)
imax = [];
imin = [];
bad = (d==0);
dd = diff([0 bad 0]);
debs = find(dd == 1);
fins = find(dd == -1);
if debs(1) == 1
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
