关于tfrrmsc问题研究
function = tfrrmsc(x,t,N,f0T,trace,K);%TFRRMSC Reassigned Morlet Scalogram time-frequency distribution.
% = TFRRMSC(X,T,N,F0T,TRACE)
% computes the Morlet scalogram and its reassigned version.
%
% X : analysed signal
% T : the time instant(s) (default : 1:length(X))
% N : number of frequency bins (default : length(X))
% F0T : time-bandwidth product of the mother wavelet
% (default : 2.5))
% TRACE : if nonzero, the progression of the algorithm is shown
% (default : 0).
% TFR,: time-frequency representation and its reassigned
% RTFR version. When called without output arguments,
% TFRRMSC runs TFRQVIEW.
% HAT : Complex matrix of the reassignment vectors.
%
% Example :
% sig=fmlin(64,0.1,0.4); tfrrmsc(sig,1:64,64,2.1,1);
%
% See also all the time-frequency representations listed in
% the file CONTENTS (TFR*)
% F. Auger, January, April 1996.
% Copyright (c) 1996 by CNRS (France).
%
% ------------------- CONFIDENTIAL PROGRAM --------------------
% This program can not be used without the authorization of its
% author(s). For any comment or bug report, please send e-mail to
% f.auger@ieee.org
if (nargin == 0), %?nargin=0?,????'At least 1 parameter required'
error('At least 1 parameter required');
end;
= size(x); a
if (nargin == 1),
t=1:xrow; N=xrow; f0T=2.5; trace=0; K= 0.001;
elseif (nargin == 2),
N=xrow; f0T=2.5; trace=0; K= 0.001;
elseif (nargin == 3),
f0T=2.5; trace=0; K=0.001;
elseif (nargin == 4),
trace = 0; K=0.001;
elseif (nargin == 5),
K= 0.001;
end;
if (N<0), error('N must be greater than zero'); end;
= size(t);
if (xcol~=1),
error('X must have only one column');
elseif (trow~=1),
error('T must only have one row');
elseif (f0T<=0),
error('F0T must be positive');
elseif (length(f0T)~=1),
error('F0T must be a scalar');
end;
if (tcol==1),
Dt=1;
else
Deltat=t(2:tcol)-t(1:tcol-1);
Mini=min(Deltat); Maxi=max(Deltat);
if (Mini~=Maxi),
error('The time instants must be regularly sampled.');
else
Dt=Mini;
end;
clear Deltat Mini Maxi;
end;
tfr= zeros(N,tcol); tf2= zeros(N,tcol);
if trace, disp('Morlet Scalogram'); end;
M=ceil(f0T*N*sqrt(2.0*log(1/K)));
tau = 0:M+round(N/2);
pi2 = 2.0*pi;
hstar = exp(-(tau/(N*f0T)).^2 /2.0) .* exp(-j*pi2*tau/N);
Thstar = tau.*hstar;
for m=1:round(N/2)-1
if trace, disprog(m,N/2,10); end;
factor=sqrt(m/(f0T*N));
for icol=1:tcol,
ti= t(icol);
tauneg=1:min();
taupos=0:min();
% positive frequencies
tfr(1+m,icol)= hstar(1+taupos*m)*x(ti+taupos);
tf2(1+m,icol)=Thstar(1+taupos*m)*x(ti+taupos);
if length(tauneg) > 0,
tfr(1+m,icol)=tfr(1+m,icol) + conj( hstar(1+tauneg*m))*x(ti-tauneg);
tf2(1+m,icol)=tf2(1+m,icol) - conj(Thstar(1+tauneg*m))*x(ti-tauneg);
end;
% negative frequencies
tfr(N+1-m,icol)=conj( hstar(1+taupos*m))*x(ti+taupos);
tf2(N+1-m,icol)=conj(Thstar(1+taupos*m))*x(ti+taupos);
if length(tauneg) > 0,
tfr(N+1-m,icol)=tfr(N+1-m,icol) +hstar(1+tauneg*m)*x(ti-tauneg);
tf2(N+1-m,icol)=tf2(N+1-m,icol) - Thstar(1+tauneg*m)*x(ti-tauneg);
end;
end;
tfr(1+m,:)=factor*tfr(1+m,:); tf2(1+m,:)=factor*tf2(1+m,:)/m;
tfr(N+1-m,:)=factor*tfr(N+1-m,:); tf2(N+1-m,:)=factor*tf2(N+1-m,:)/m;
end;
m=round(N/2);
factor=sqrt(m/(f0T*N));
if trace, disprog(m,N/2,10); end
for icol=1:tcol,
ti= t(icol);
tauneg=1:min();
taupos=0:min();
tfr(1+m,icol)= hstar(1+taupos*m)*x(ti+taupos);
tf2(1+m,icol)=Thstar(1+taupos*m)*x(ti+taupos);
if length(tauneg) > 0,
tfr(1+m,icol)=tfr(1+m,icol) + conj( hstar(1+tauneg*m))*x(ti-tauneg);
tf2(1+m,icol)=tf2(1+m,icol) - conj(Thstar(1+tauneg*m))*x(ti-tauneg);
end;
end;
tfr(1+m,:)=factor*tfr(1+m,:);
tf2(1+m,:)=factor*tf2(1+m,:)/m;
avoid_warn=find(tfr~=0.0);
tf2(avoid_warn)=tf2(avoid_warn)./tfr(avoid_warn);
tfr=abs(tfr).^2;
if trace, disp('reassignment :'); end;
rtfr= zeros(N,tcol);
Ex=mean(abs(x(min(t):max(t))).^2); Threshold=1.0e-6*Ex;
factor=2.0*pi*N*f0T*f0T;
for icol=1:tcol,
if trace, disprog(icol,tcol,10); end;
for jcol=1:N,
if tfr(jcol,icol)>Threshold,
icolhat= icol + round(real(tf2(jcol,icol)/Dt));
icolhat=min(max(icolhat,1),tcol);
m=rem(jcol+round(N/2)-2,N)-round(N/2)+1;
jcolhat= jcol + round(imag(m*m*tf2(jcol,icol)/factor));
jcolhat=rem(rem(jcolhat-1,N)+N,N)+1;
rtfr(jcolhat,icolhat)= rtfr(jcolhat,icolhat)+tfr(jcol,icol);
tf2(jcol,icol)= jcolhat + 1j * icolhat;
else
tf2(jcol,icol)=(1+j)*inf;
rtfr(jcol,icol)=rtfr(jcol,icol)+tfr(jcol,icol);
end;
end;
end;
if (nargout==0),
continu=1;
while (continu==1),
choice=menu ('Choose the representation:',...
'stop',...
'Morlet scalogram',...
'reassigned Morlet scalogram');
if (choice==1), continu=0;
elseif (choice==2),
tfrqview(tfr,x,t,'tfrmsc',f0T);
elseif (choice==3),
tfrqview(rtfr,x,t,'tfrrmsc',f0T);
end;
end;
elseif (nargout>2),
hat=tf2;
end;
谁能解析下这段程序啊、还有tfr与rtfr的区别啊
页:
[1]