insttf_abs.m

Matlab时频分析工具箱,希望能对大家有所帮助啊

function [itime,ifreq]=insttf_abs(s,g,a,difforder)
%INSTFREQ   Instantaneous frequency computed from absolute values.
%   Usage:  [itime,ifreq] = insttf_abs(s,g,a);
%           [itime,ifreq] = insttf_abs(s,g,a,order);
%
%   [itime,ifreq]=INSTFTF_ABS(s,g,a) computes the instantanous time and
%   frequency from the spectrogram s. The spectrogram must have been
%   computed using the window g and time-shift a.
%
%   [itime,ifreq]=INSTFTF_ABS(s,g,a,difforder) uses a centered finite
%   diffence scheme of order difforder to perform the needed numerical
%   differentaion. Default is to use a 4th order scheme.
%
%   Currently this method only works if the window g is 'gauss' or a cell
%   array {'gauss',tfr} where tfr is the time-to-frequency ratio of the
%   Gaussian.
%
%   See the help text on INSTTF_DGT for a description of the output itime
%   and ifreq.
%
%   If the original signal is known, the instantaneous time and frequency
%   can be more accuratly computed by INSTTF_DGT
% 
%   SEE ALSO: INSTTF_DGT, INSTTF_PHASE, REASSIGN
% 
%   REFERENCES:
%     E. Chassande-Mottin, I. Daubechies, F. Auger, and P. Flandrin.
%     Differential reassignment. Signal Processing Letters, IEEE,
%     4(10):293-294, 1997.

% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.
  
% Author: Peter L. Soendergaard, 2008.

error(nargchk(3,4,nargin));

if ~isreal(s)
  error('First input argument must be real-valued.');
end;

if ~(all(s(:)>=0))
  error('First input argument must be positive or zero.');
end;
  
if nargin==3
  difforder=4;
end;

M=size(s,1);
N=size(s,2);
W=size(s,3);

L=N*a;
tfr=1;

if isnumeric(g)
  error('Only Gaussian windows supported so far.');
end;

[g,info]=comp_window(g,a,M,L,0);

if ~info.gauss
  error(['The window must be a Gaussian window (specified as a string or ' ...
         'as a cell arrray).']);
end;

L=N*a;
b=L/M;

% We must avoid taking the log of zero.
% Therefore we add the smallest possible
% number
logs=log(s+realmin);
  
itime=pderiv(logs,2,difforder)/(2*pi)*info.tfr;
ifreq=pderiv(logs,1,difforder)/(2*pi*info.tfr);