tfrunter.m

时间序列的时频分析工具箱

function [tfr,t,f]=tfrunter(X,time,form,fmin,fmax,N,trace); 
%TFRUNTER Unterberger time-frequency distribution.
%	[TFR,T,F]=TFRUNTER(X,T,FORM,FMIN,FMAX,N,TRACE) generates the
%	auto- or cross-Unterberger distribution (active or passive form).
%
%	X : signal (in time) to be analyzed. If X=[X1 X2], TFRUNTER 
%	   computes the cross-Unterberger distribution (Nx=length(X)).
%	T : time instant(s) on which the TFR is evaluated  
%	    					(default : 1:Nx).
%	FORM : 'A' for active, 'P' for passive Unterberger distribution.
%						(default : 'A'). 
%	FMIN,FMAX : respectively lower and upper frequency bounds of 
%	   the analyzed signal. These parameters fix the equivalent 
%	   frequency bandwidth (expressed in Hz). When unspecified, you
%	   have to enter them at the command line from the plot of the
%	   spectrum. FMIN and FMAX must be >0 and <=0.5. 
%	N : number of analyzed voices (default : automatically determined).
%	TRACE : if nonzero, the progression of the algorithm is shown
%	                                	(default : 0).
%	TFR : time-frequency matrix containing the coefficients of the 
%	   decomposition (abscissa correspond to uniformly sampled
%	   time, and ordinates correspond to a geometrically sampled 
%	   frequency). First row of TFR corresponds to the lowest 
%	   frequency. When called without output arguments, TFRUNTER
%	   runs TFRQVIEW.
%	F : vector of normalized frequencies (geometrically sampled 
%	   from FMIN to FMAX).
% 
%	Example : 
%	 sig=altes(64,0.1,0.45); tfrunter(sig);
%
%	See also all the time-frequency representations listed in
%	the file CONTENTS (TFR*)

%	P. Goncalves, November 95 - O. Lemoine, June 1996.
%	Copyright (c) 1995 Rice University - CNRS (France) 1996.
%
%	------------------- 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 
%	[e-mail] 

global ratio_f

if (nargin == 0),
 error('At least one parameter required');
end;

[xrow,xcol] = size(X);
if nargin<=6, trace=0; end

if (nargin == 1),
 time=1:xrow; form='A'; 
elseif nargin==2,
 form='A'; 
elseif (nargin==4),
 disp('FMIN will not be taken into account. Determine it with FMAX');
 disp('     from the following plot of the spectrum.'); 
elseif nargin==5,
 N=[];
end;

form=upper(form);

[trow,tcol] = size(time);
if (xcol==0)|(xcol>2),
 error('X must have one or two columns');
elseif (trow~=1),
 error('T must only have one row'); 
elseif (form~='A' & form~='P'),
 error('Unknown value for FORM'); 
end; 

Mt = length(X); 

if trace, disp('Unterberger distribution'); end;

if xcol==1,
 X1=X;
 X2=X; 
else
 X1=X(:,1);
 X2=X(:,2);
end
s1 = real(X1);
s2 = real(X2);
M  = (Mt+rem(Mt,2))/2;

t = (1:Mt)-M-1;
T = xrow;

if nargin<=3,				        % fmin,fmax,N unspecified
 STF1 = fft(fftshift(s1(min(time):max(time)))); Nstf=length(STF1);
 sp1 = (abs(STF1(1:Nstf/2))).^2; Maxsp1=max(sp1);
 STF2 = fft(fftshift(s2(min(time):max(time)))); 
 sp2 = (abs(STF2(1:Nstf/2))).^2; Maxsp2=max(sp2);
 f = linspace(0,0.5,Nstf/2+1) ; f=f(1:Nstf/2);
 plot(f,sp1) ; grid; hold on ; plot(f,sp2) ; hold off
 xlabel('Normalized frequency');
 title('Analyzed signal energy spectrum');
 axis([0 1/2 0 1.2*max(Maxsp1,Maxsp2)]) ; 
 indmin=min(find(sp1>Maxsp1/100));
 indmax=max(find(sp1>Maxsp1/100));
 fmindflt=max([0.01 0.05*fix(f(indmin)/0.05)]);
 fmaxdflt=0.05*ceil(f(indmax)/0.05);
 txtmin=['Lower frequency bound [',num2str(fmindflt),'] : '];
 txtmax=['Upper frequency bound [',num2str(fmaxdflt),'] : '];
 fmin = input(txtmin); fmax = input(txtmax);
 if isempty(fmin), fmin=fmindflt; end
 if isempty(fmax), fmax=fmaxdflt; end
end

if fmin >= fmax
 error('FMAX must be greater or equal to FMIN');
elseif fmin<=0.0 | fmin>0.5,
 error('FMIN must be > 0 and <= 0.5');
elseif fmax<=0.0 | fmax>0.5,
 error('FMAX must be > 0 and <= 0.5');
end

B = fmax-fmin ; 
R = B/((fmin+fmax)/2) ; 
ratio_f = fmax/fmin ;    
umax = fzero('umaxunte',0); 
Teq = M/(fmax*umax);  
if Teq<2*M,
 M0 = round((2*M^2)/Teq-M)+1;
 M1 = M+M0;
 T  = 2*M1-1;
elseif Teq>=2*M
 M0 = 0;
 M1 = M;
end;

Nq= ceil((B*T*(1+2/R)*log((1+R/2)/(1-R/2)))/2);
Nmin = Nq-rem(Nq,2);
Ndflt = 2^nextpow2(Nmin);
if nargin<=3,
 Ntxt=['Number of frequency samples (>=',num2str(Nmin),') [',num2str(Ndflt),'] : '];
 N = input(Ntxt);
end
if ~isempty(N),
 if (N<Nmin),
  dispstr=['Warning : the number of analyzed voices (N) should be > ',num2str(Nmin)];
  disp(dispstr);
 end
else
 N=Ndflt; 
end

fmin_s = num2str(fmin) ; fmax_s = num2str(fmax) ; N_s = num2str(N) ;
if trace,
 disp(['frequency runs from ',fmin_s,' to ',fmax_s,' with ',N_s,' points']);
end


% Geometric sampling of the analyzed spectrum
k = 1:N;
q = (fmax/fmin)^(1/(N-1));
t = (1:Mt)-M-1;
geo_f  = fmin*(exp((k-1).*log(q)));
tfmatx = zeros(Mt,N);
tfmatx = exp(-2*i*t'*geo_f*pi);
S1 = s1'*tfmatx; 
S2 = s2'*tfmatx; 
clear tfmatx
S1(N+1:2*N) = zeros(1,N);
S2(N+1:2*N) = zeros(1,N);


% Mellin transform computation of the analyzed signal
p = 0:(2*N-1);
Mellin1 = fftshift(ifft(S1));
Mellin2 = fftshift(ifft(S2));
umin = -umax;
du = abs(umax-umin)/(2*M1);
u(1:2*M1) = umin:du:umax-du;
u(M1+1) = 0;
beta = (p/N-1)./(2*log(q));


% Computation of the Lambda(+/- u) dilations/compressions 
% of the analyzed signal
waf = zeros(2*M1,N);
for n = 1:2*M1,
 if trace, disprog(n,4*M1,10); end
 MX1 = exp(-(2*i*pi*beta+0.5)*log(sqrt(1+(u(n)/2)^2)-u(n)/2)).*Mellin1 ;
 MX2 = exp(-(2*i*pi*beta+0.5)*log(sqrt(1+(u(n)/2)^2)+u(n)/2)).*Mellin2 ;
 FX1 = fft(fftshift(MX1)) ;
 FX1 = FX1(1:N) ;
 FX2 = fft(fftshift(MX2)) ;
 FX2 = FX2(1:N) ;
 if form=='A',
  waf(n,:) = FX1.*conj(FX2); 
 else
  waf(n,:) = FX1.*conj(FX2)*(1/sqrt(1+(u(n)/2)^2)); 
 end
end
waf  = [waf(M1+1:2*M1,:) ; waf(1:M1,:)].*geo_f(ones(2*M1,1),:);
tffr = ifft(waf);  
tffr = real(rot90([tffr(M1+1:2*M1,:) ; tffr(1:M1,:)],-1));


% Conversion from [t.f,f] to [t,f] using a 1-D interpolation
tfr  = zeros(N,tcol);
Ts2  = (Mt-1)/2 ;
gamma = linspace(-geo_f(N)*Ts2,geo_f(N)*Ts2,2*M1) ;
alpha = (0.6*N-1)/0.4;
for n = 1:N,
 if trace, disprog(n+alpha,N+alpha,10); end
 ind = find(gamma>=-geo_f(n)*Ts2 & gamma<=geo_f(n)*Ts2);
 x = gamma(ind);
 y = tffr(n,ind);
 xi = (time-Ts2-1)*geo_f(n);
 vnum=version;
 if vnum(1)=='5',
  tfr(n,:) = interp1(x,y,xi,'spline');
 else
  tfr(n,:) = interp1(x,y,xi,'spline')';
 end;
 clear vnum
end 


t = time;
f = geo_f(1:N)';

% Normalization
SP1 = fft(hilbert(s1)); 
SP2 = fft(hilbert(s2)); 
indmin = 1+round(fmin*(xrow-2));
indmax = 1+round(fmax*(xrow-2));
SP1ana=SP1(indmin:indmax);
SP2ana=SP2(indmin:indmax);
tfr=tfr*(SP1ana'*SP2ana)/integ2d(tfr,t,f)/N;

clear ratio_f

if (nargout<=1),
 tfrqview(tfr,hilbert(real(X1)),t,'TFRUNTER',N,f);
end;