nfame.m

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

function [nf, na, ex, E, N] = nfame(data,dt)
%
%    [nf,na,ex,E,N]=nfame(data,dt) :
%
%    Function to calculate the modified Hilbert frequency and amplitude
%    of data(n,k), where n specifies the length of time series, and k is 
%    the IMF component number.
%    Based on normalized data with spline fitting of maxima with the ends fixed.
%    Final frequency and amplitude values are smoothed with a five point median filter.
%    Modified Hilbert transform function 'hilbtm' with the end effect
%    eliminated by extension of data is used to perform a Hilbert transform.
%
%    Input-
%	data	- 2-D matrix data(n,k), that specifies one IMF component,
%		  where n is the number of points, and 
%		  k is the IMF component number
%	dt	- sampling period in seconds
%    Output-
%	nf	- 2-D matrix nf(n,k) that specifies the Hilbert frequency in Hz
%	na	- 2-D matrix na(n,k) that specifies the Hilbert Amplitude
%	ex	- 2-D matrix ex(n,k) that specifies the Splined envelope
%	E	- 2-D matrix e(n,k) that specifies the Normalized input data
%	N	- 2-D matrix N(n,k) that specifies the Nonlinearity Index
%
%    Norden Huang  (NASA GSFC)	June 2, 2002 Initial
%    J.Marshak (NASA GSFC)    November 8, 2003 Edited the description
%
%    Notes-
%    Non MATLAB Library routines used in the function are:
%	emax.m
%	hilbtm.m
%
%    Temporary remarks-
%    1) Size of tx and mx is the same.
%    2) The code is the same as nfam5, except for the parameter naming and modification date.
%    3) Added "and amplitude".
%    4) Works only for one IMF component.

%----- Get dimensions
[n,m] = size(data); 
te=1:n; te=te';

%----- Extract the set of max points and their coordinates
[mx, tx]=emax(data);

%----- Fix the ends to prevent wide swaying in spline 
%----- by adding the te(1) anf te(n) and mx(1) and mx(p) 
%----- to the first and last tx and mx
[p,q]=size(mx);
[u,v]=size(tx);
tx=[te(1);tx;te(n)];
mx=[mx(1);mx;mx(p)];
ex=spline(tx,mx,te);

%----- Normalize data by splined envelope
E=data./ex;

%----- Apply modified Hilbert transform
nf=diff(unwrap(angle(hilbtm(E))))./(2*pi*dt);
nf=[nf;nf(n-1)];
na=abs(hilbtm(E));

%----- Apply median for every 5 points, 5/4/02
for i=3:n-2
    nf(i,:)=median(nf(i-2:i+2,:));
    na(i,:)=median(na(i-2:i+2,:));
end
    
%----- Compute the Nonlinear Index
N=(na-1).^2;

%----- Re-constitute na
%  na=na.*ex;

%----- Plot the results
plot(te,data,te,ex,te,E,te,nf,te,na,te,N, 'LineWidth', 1.5);
legend('Data','ex','E','nf','na','N');