function h=hilbtm(x)
% hilbtm(x): Improved Hilbert transform on x(n,k)
%        k: the sifted components; n: length of the time series.

%	Z. SHEN    Jan. 1996   Initial
%  D. XIANG   04-04-2002  Modified At the Johns Hopkins University
%  Karin Blank 4/18/03 Modified

[n,k]=size(x);
for j=1:k
    zz= [];

    %Modified code to correctly detect extrema -KB
    n_mx = 1; n_mn = 1;
    flagmax=0; flagmin=0;
    for i=2:n-1
        %finds local maximum
        if ((x(i,j) > x(i-1,j)) & (flagmax == 0))
            flagmax = 1;
            Xmax = i;
        end    
        if((x(i,j) < x(i+1,j)) & (flagmax == 1))
            flagmax = 0;
        end
        if((x(i,j) > x(i+1,j)) & (flagmax == 1))
            flagmax = 0;
            n_mx = Xmax;
        end
        
        %finds local minimum
        if((x(i,j) < x(i-1,j)) & (flagmin == 0))
            flagmin = 1;
            Xmin = i;
        end
        if((x(i,j) > x(i+1)) & (flagmin == 1))
            flagmin = 0;
        end
        if((x(i) < x(i+1)) & (flagmin == 1))
            flagmin = 0;
            
            n_mn = Xmin;
        end 
        
        if(x(n_mn) < 0 & x(n_mx) > 0) & (n_mn > 1 & n_mx > 1) %if 2 extrema found and they are across the zero crossing
            break;
        end
    end
    %end extrema detection
    
    if n_mn<n_mx,
        first = n_mn;
        second = n_mx;
    else
        first = n_mx;
        second = n_mn;
    end
    
    zz=x(first:second,j); %get all y values between min and max extrema
    mm=length(zz); %length of zz
    ia=fix(first/(2*(mm-1)))+1; %calculate the number of times to copy artificial wave
    iaa=max(ia,2); %DX, modify constant from 1 to 2. %do it at least twice
    zz1=flipud(zz); %zz1 is zz (all the y values between the min and max extrema) reversed
    
    if n_mn<n_mx
        x1=zz1; %copy those values to x1
    else
        x1=[zz;zz1(2:mm)];
    end
    for jj=1:iaa,
        x1=[x1;zz(2:mm);zz1(2:mm)];  % copy the artificial wave a couple times
    end
    x1(:,j)=x1;
    %DX, find the first zero-cross and the slop
    r=length(x1);
    for kk=1:r-1
        if((x1(kk,j)*x1(kk+1,j))<=0)
            if(abs(x1(kk,j)) > abs(x1(kk+1,j)))
                n0=kk; %found zero crossing
            else
                n0=kk+1; %found zero crossing
            end
            s0=x1(kk+1,j)-x1(kk,j); %calculate slope
            break;
        end
    end
    x1=x1(n0:r,j); %copy all artificial wave values from zero crossing to end 
    
    x1=[x1;x(first+1:n,j)]; %attach actual y values from the minimum extrema to the end of the dataset
    sz=length(x1); %length of x1
    np1=sz-n; %length of artificial wave portion of x1

 

    %---------------------Treat the tail----------------------------

    %Modified code to correctly detect extrema -KB
    n_mx = 1; n_mn = 1;
    flagmax=0; flagmin=0;
    x_flip = fliplr(x1');
    for i=2:sz-1
        %finds local maximum
        if ((x_flip(i) > x_flip(i-1)) & (flagmax == 0))
            flagmax = 1;
            Xmax = sz-i;
        end    
        if((x_flip(i) < x_flip(i+1)) & (flagmax == 1))
            flagmax = 0;
        end
        if((x_flip(i) > x_flip(i+1)) & (flagmax == 1))
            flagmax = 0;
            n_mx = Xmax;
        end
        
        %finds local minimum
        if((x_flip(i) < x_flip(i-1)) & (flagmin == 0))
            flagmin = 1;
            Xmin = sz-i;
        end
        if((x_flip(i) > x_flip(i+1)) & (flagmin == 1))
            flagmin = 0;
        end
        if((x_flip(i) < x_flip(i+1)) & (flagmin == 1))
            flagmin = 0;
            
            n_mn = Xmin;
        end 
        
        if(x_flip(sz-n_mn) < 0 & x_flip(sz-n_mx) > 0) & (n_mn > 1 & n_mx > 1) %if 2 extrema found and they are across the zero crossing
            break;
        end
    end
    clear x_flip;
    %end extrema detection

    if n_mn<n_mx
        first = n_mn;
        second = n_mx;
    else
        first = n_mx;
        second = n_mn;
    end
    zz=x1(first:second); %use all y values between xmin and xmax
    mm=length(zz); 
    zz1=flipud(zz);	%reverse order y values
    ia=fix((sz-second)/(2*(mm-1)))+1; %calculate number of copies of artificial wave to make
    iaa=max(ia,2); %DX, modify constant from 1 to 2. %at least 2
    
    if(n_mn< n_mx)       
        x1=[x1(1:second);zz1(2:mm)]; %attach initial wave segment
        for jj=1:iaa
            x1=[x1;zz(2:mm);zz1(2:mm)]; %copies artificial waves and appends them to x1
        end
        x1=[x1;zz(2:mm-1)]; %attaches last wave segment
    else
        x1=x1(1:second);
        for jj=1:iaa
            x1=[x1;zz1(2:mm);zz(2:mm)]; %copies artificial waves and appends them to x1
        end
        
        x1=[x1;zz1(2:mm-1)]; %attaches last wave segment
    end
    x1(:,j)=x1;


    
    %DX, find the first zero-cross and the slop
    r = length(x1);
    for k=1:r-1
        kk=r-k+1;
        if((x1(kk,j)*x1(kk-1,j))<=0)
            if(abs(x1(kk,j)) > abs(x1(kk-1,j)))
                n0=kk;
            else
                n0=kk-1;
            end
            
            if(((x1(kk,j)-x1(kk-1,j))*s0)>0)   %same sign as s0
                break;
            end
        end
    end
    x1=x1(1:n0,j);


	x1=hilbert(x1); %calculate hilbert transform
    if ((np1+n-1)<=n0)
	    x(:,j)=x1(np1:np1+n-1); %truncate artificial elements of dataset
    else
        np2=n0-np1;
        x(1:np2,j)=x1(np1:n0-1);
    end
end
h=x;
clear x