fa.m

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

function [freq,amp] = fa(c,ifmethod,normmethod,nfilter)
% COMPONENTS/FA finds the frequency and amplitude of components
%
% Usage:
%  [freq,amp] = fa(c,ifmethod,normmethod,nfilter)
% All input parameters after 'c' are optional.
%
% Input:
%  c          - the Components data structure containing the input data
%  ifmethod   - Handle to the function to call to get the freqency and
%               amplitude of a component. See notes.
%  normmethod - Handle to a function to call to normalize a component, or
%               'none'. See notes.
%  nfilter    - Number of points to smooth the data with a median filter,
%               or 0 to skip smoothing (recommended)
% Output: (all outputs optional)
%  freq       - A Components object holding the instantaneous frequency for
%               each component in 'c' except the trend
%  amp        - A Components object holding the instantaneous amplitude for
%               each component in 'c' except the trend
%
% Notes:
%    MATLAB function handles have the syntax '@function_name'. So to use the
%  default ifmethod, the Hilbert transform, pass '@fah' (without quotes) as
%  the ifmethod parameter.
%
%    The ifmethod function passed must accept the following prototype:
%  [f,a] = faX (data, dt)
%  where 'data' is the vector of the current component, dt is delta-t (the
%  increment in time per sample), f is the instantaneous frequency of the
%  component, and a is the instantaneous amplitude. f and a must have the
%  same length as 'data'.
%
%    The following ifmethod functions are distributed with the HHT library:
%   - fah: Hilbert transform (default) (normalization helps; see below)
%   - faz: Huang's Generalized Zero Crossing algorithm (do not normalize)
%   - faacos: Arccosine-based instantaneous frequency (MUST be normalized)
%   - desa*: Teager Energy Operator (normalization effects unstudied)
%
%    The normmethod function passed must accept the following prototype:
%  [nd,na] = Xnormalize (data)
%  where 'data' is the vector of the current component, nd is the
%  normalized data, and na is the amplitude by which the data was
%  normalized, i.e., nd = data ./ na.
%    The following normmethod functions are distributed with the HHT
%  library:
%   - splinenormalize: estimates data envelope by data extrema,
%     interpolated by a cubic spline (default when ifmethod is @fah or @faacos)
%   - hilbertnormalize: uses Hilbert transform to get envelope (more
%     sensitive to intra-wave ripples)
%   - blocknormalize: normalize by assuming constand amplitude between zero
%     crossings (not recommended)
%    As a special case, pass 'none' as normmethod to prevent any
%  normalization, even the default.
%
% Examples: ('c' holds the components of interest)
%
% * Get the frequency and amplitude using the default method (normalized Hilbert):
%  >> [freq,amp] = fa(c);
%
% * Use Huang's Generalized Zero Crossing method:
%  >> [freq,amp = fa(c,@faz);

% Kenneth C. Arnold (for NASA GSFC), 2004-08-06

if nargin<2; ifmethod=[]; end
if nargin<3; normmethod=[]; end
if nargin<4; nfilter=[]; end

if isempty(ifmethod); ifmethod = @fah; end
if isempty(normmethod) & (isequal(ifmethod, @fah) | isequal(ifmethod, @faacos))
    normmethod = @splinenormalize;
end
if isequal(normmethod,'none'); normmethod=[]; end
if isempty(nfilter); nfilter=0; end

freq=c;
amp=c;

% cut off the space for the residual component
freq.d(end) = [];
amp.d(end) = [];

nc=get(c,'nc');
for i=1:nc-1 % all but the residual component
    %--- Normalize data, if requested.
    if ~isempty(normmethod)
        [nd,na]=feval(normmethod,c.d(i).c);
    else
        nd = c.d(i).c;
    end
    
    %--- Calculate frequency and amplitude
    [freq.d(i).c,amp.d(i).c] = feval(ifmethod, nd, get(c,'dt',i));

    %--- Replace amplitude from ifmethod with normalization amplitude
    if ~isempty(normmethod)
        amp.d(i).c = na;
    end
    
    %--- Filter data, if requested.
    if nfilter > 0
        freq.d(i).c = medianfilter(freq.d(i).c,nfilter);
    end
end