splinenormalize.m

这是一个关于hht变换很有用的工具箱

function [data,a]=splinenormalize(data)
 
% The function SPLINENORMALIZE normalizes with splines,
% the data(n,m) where n specifies the number of time points, 
% and m is the number of IMF components.
%
% Calling sequence-
% [data,a]=splinenormalize(data)
%
% Input-
%	data	- 2-D matrix data(n,m) that specifies the IMF components 
% Output-
%	data	- normalized data
%	a	    - splined envelope
%
% Used by-
%	FA
 
% K. Arnold (for NASA GSFC)	    Jan. 28 2004 Modified
%				(fixed a bug where a "random" point 
%				could get thrown somewhere in the middle
%				of the dataset and the endpoints 
%				were not controlled).
% K. Arnold (for NASA GSFC)     Aug. 4 2004 Modified
%  (made more accurate by using algorithm from extrema.m to find the extrema)
% K. Arnold (for NASA GSFC)     Aug. 5 2004 Modified
%  (speed up extrema loop by avoiding appending, which is expensive in MATLAB)

%----- Get the dimension
[n,m] = size(data);

%----- Initialize and flip data if needed 
flipped=0;
if (n<m)
    flipped=1;
    data=data';
    [n,m]=size(data);
end

te=(1:n)';

%----- Process each IMF component
for c=1:m
    %----- Extract the data extrema (see 'extrema.m')
    x=zeros(1,n-2);
    y=x;
    nExtrema = 0;
    wasSmaller = data(2,c) > data(1,c);
    for i=2:n-1
        isSmaller = data(i+1,c) > data(i,c);
        if wasSmaller ~= isSmaller & data(i+1,c) ~= data(i,c)
            diff = data(i+1,c)-data(i-1,c);
            den = (2*data(i,c)-data(i-1,c)-data(i+1,c));
            nExtrema = nExtrema+1;
            x(nExtrema) = i + diff / (2*den);
            y(nExtrema) = abs(data(i,c) + (diff*diff) / (8*den));
            wasSmaller = isSmaller;
        end 
    end
    
    %----- Fix the end to prevent wide swaying in spline 
    %----- by assigning the te(1) and te(n) 
    %----- the same values as the first and last tx and mx.
    if nExtrema > 1
        x=[1 x(1:nExtrema) n];
        %----- Fix the ends at the same as the next point
        y=[y(1) y(1:nExtrema) y(nExtrema)];
        a(:,c)=spline(x,y,te);
        %----- %  Normalize the data by splined envelope
        data(:,c)=data(:,c)./a(:,c);
    else
        %----- Leave data unchanged
        a(:,c)=ones(n,1);
    end
end

%----- Flip data back if needed
if (flipped)
    data=data';
    a=a';
end