ifp.m

用于信号处理里面的时频分析

function xp1=ifp(x,dt,epsi)
n=length(x);
x1=hilbert(x);
xx=real(x1);
yy=imag(x1);
gmax=0.0;
for jj=1:n
    aaa=xx(jj)*xx(jj)+yy(jj)*yy(jj);
	if aaa>gmax
        gmax=aaa;
    end
end
for ii=1:n
    if ii>1 & ii<n-1
    ss1=xx(ii+1)-xx(ii-1);
	ss1=ss1/dt/2.;
	ss2=yy(ii+1)-yy(ii-1);
	ss2=ss2/dt/2.;
	aaa=xx(ii)*xx(ii)+yy(ii)*yy(ii);
	xp1(ii)=xx(ii)*ss2-yy(ii)*ss1;
	xp1(ii)=xp1(ii)/(aaa+gmax*epsi)/2/pi;
    if xp1(ii)<0.0
        xp1(ii)=0.0;
    end
else
    xp1(ii)=0.0;
end
end