maxp.m

一种新的时频分析方法的matlab源程序。

function[x,y,xx,yy]=maxp(nlook,data,nshift,periodk,periodk1)

% The function MAXP extracts maximum points and their coordinates for the data with well distinguishable extrema. 
% The data do not have to be preprocessed by EMD.
% The procedure includes finding the defined number of descending global 
% maximum points, analyzing there period and taking the maximum of them.
% The maximum data values are searched within each so defined running
% periods, saved on screen and to the MATLAB format file, and plotted.
%
% First two arguments are required, others are optional, but if used
% should appear in the order defined in the parameter list.
%
% Example
% 	[x,y,xx,yy]=maxp(100,AS1A);
%	The missing parameters are assigned by default.
%
% Calling sequence-
% [x,y,xx,yy]=maxp(nlook,data[,nshift][,periodk][,periodk1])
%
% Input-
%	data	- 2-D matrix of data that specifies the data values and
%		    their coordinates
%	nlook	- integer number that specifies the number of initial 
%		    global maximums that will be found in descending order
%	nshift	- integer number that specifies the number of points
%		    that defines how much each period beginning 
%		    is shifted from the previous extrema
%	periodk	- the koefficent that specifies the extent of the period
%	periodk1- the koefficent that specifies the extent of the first
%		    period
%
% Output-
%	x	    - vector that specifies the coordinates of data
%	y	    - vector that specifies the values of data
%	xx	    - vector that specifies the coordinates of data maximums
%	yy	    - vector that specifies the maximum values of data

% J.Marshak (NASA GSFC)	Dec.22, 2003 Initial

%----- Define default parameters
if nargin < 1
   error('maxp: first 2 parameters must be passed');
end
if nargin < 5
    periodk1=1.1;
end
if nargin < 4
    periodk=1.5;
end
if nargin < 3
    nshift=30;
end
%if nargin < 2
%    load C:\HHT\lena\Data\AB1A.TXT;
%    data=AB1A;
%end


%----- Get the X-coordinates and their Y-values
x=data(:,1);
y=data(:,2);
m=size(x,1);

%----- Find first nlook global maximum points
zm=y;
gdmxi=1;
for ilook=1:nlook
    gdmxy=-999;
    for i=1:m
        if (zm(i) > gdmxy & zm(i) ~= 999.)
           gdmxi=i;
           gdmxx=x(i);
           gdmxy=y(i);
        end
    end
    zm(gdmxi)=999.;
end

%----- Define the period where to look for maximum
%----- and print a few
zx(1:nlook)=0;
zy(1:nlook)=0;
zp(1:nlook)=0;
zi(1:nlook)=0;
zp(1)=0.;
j=0;

fprintf(1,'       Order                Index Period    X      Y \n')
for i=1:m
    if (zm(i) == 999.)
       j=j+1;
       zy(j)=y(i);
       zx(j)=x(i);
       zi(j)=i;
       if ( j > 1 ) 
           zp(j)=zx(j)-zx(j-1);
       end
       if ( j < 200 )
           line=[zi(j);zp(j);zx(j);zy(j)];
           fprintf(1,'%8d LINE OF TABLE %8d,%6.2f,%6.2f,%6.2f\n',j,line)
       end
    end
end
maxPeriod=max(zp);
mxp=maxPeriod*periodk1;
mxpr=maxPeriod*periodk;
fprintf(1,'\nMAX Period: %6.2f Extended Period: %6.2f\n\n',maxPeriod,mxpr);

%----- Find maximum points and their coordinates within each defined period
tmpx(1:m)=0.;
tmpy(1:m)=0.;
ib=1;

for jmx=1:m
    jdim=jmx;
    zmxy=-999.;
    for i=ib:m-1
        ie=i+1;
        period=x(ie)-x(ib);
        if( ib == 1 ) 
           mxzp=mxp;
        else 
           mxzp=mxpr;
        end
        if ( period <= mxzp )
           %----- find max within the extended period
           if ( y(i) > zmxy )
              zmxy=y(i);
              zmxx=x(i);
              zmxi=i;
           end
        else
           %----- leave the loop 
           break;
        end
     end

     %----- save the maximum
     tmpx(jmx)=zmxx;
     tmpy(jmx)=zmxy;

     %----- check the bounds
     if ( ie == m )
        break;
    else
        line=[jmx;zmxi;zmxx;zmxy;ib;ie;ie-ib+1];
        fprintf(1,'PERIOD:%8d,%8d,%6.2f,%6.2f,%8d,%8d,%8d\n',line);
        %----- define the beginning of next period
        ib=zmxi+nshift;
        if ((ib+1) <= m )
          continue;
        end
     end
 end

%----- check the square distance for the last period
difPrevious=(tmpx(jdim-1)-tmpx(jdim-2) - max(zp))^2;
difLast=(tmpx(jdim)-tmpx(jdim-1) - max(zp))^2;
if (difLast > 0.001 )
   fprintf(1,'\nSquare Distance for the accepted period   :%10.6f',difPrevious);
   fprintf(1,'\nSquare Distance for the disregarded period:%10.6f\n',difLast); 
   jdim=jdim-1;
end

%----- save the final output to the array and in file
clear xx yy;
xx(1:jdim)=0.;
yy(1:jdim)=0.;
xx(1:jdim)=tmpx(1:jdim);
yy(1:jdim)=tmpy(1:jdim);

save tmp.txt xx yy -ascii

%----- plot the data
plot(x,y,xx,yy,'r*','LineWidth', 1.5);