tfrbert.m

matlab时间序列分析工具箱

function  [tfr,t,f]=tfrbert(X,time,fmin,fmax,N,trace);
%TFRBERT Unitary Bertrand time-frequency distribution.
%	[TFR,T,F]=TFRBERT(X,T,FMIN,FMAX,N,TRACE) generates the
%	auto- or cross- unitary Bertrand distribution. 
%
%	X : signal (in time) to be analyzed. If X=[X1 X2], TFRBERT 
%	   computes the cross-unitary Bertrand distribution (Nx=length(X)).
%	T : time instant(s) on which the TFR is evaluated
%	   					(default : 1:Nx).
%	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
%	   distribution (x-coordinate corresponds to uniformly sampled 
%	   time, and y-coordinate corresponds to a geometrically sampled
%	   frequency). First row of TFR corresponds to the lowest 
%	   frequency. When called without output arguments, TFRBERT
%	   runs TFRQVIEW.
%	F : vector of normalized frequencies (geometrically sampled 
%	   from FMIN to FMAX).
%
%	Example :    
%	 sig=altes(64,0.1,0.45); tfrbert(sig);
%
%	See also TFRSCALO, TFRUNTER, TFRDFLA, TFRASPW.	
%	P. Goncalves, October 95 - O. Lemoine, June 1996. 
%	Copyright (c) 1995 Rice University - 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
%	[e-mail] 

global ratio_f

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

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

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

[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'); 
end; 

Mt = length(X); 

if trace, disp('Unitary Bertrand 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;
Tmin = 1; 
Tmax = Mt; 
T = Tmax-Tmin;

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('umaxbert',0); 
Teq = M/(fmax*umax);  
if Teq<Mt,
 M0 = round((2*M^2)/Teq-M)+1;
 M1 = M+M0;
 T  = 2*M1-1;
elseif Teq>=Mt
 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:N*2) = zeros(1,N);
S2(N+1:N*2) = 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 P0(t.f,f) 
waf = zeros(2*M1,N);
for n = [1:M1,M1+2:2*M1],
 if trace, disprog(n,4*M1,10); end
 MX1 = exp((-2*i*pi*beta+0.5)*log((u(n)/2)*exp(-u(n)/2)/sinh(u(n)/2))).*Mellin1;
 MX2 = exp((-2*i*pi*beta+0.5)*log((u(n)/2)*exp( u(n)/2)/sinh(u(n)/2))).*Mellin2;
 FX1 = fft(fftshift(MX1)) ;
 FX1 = FX1(1:N) ;
 FX2 = fft(fftshift(MX2)) ;
 FX2 = FX2(1:N) ;
 waf(n,:) = FX1.*conj(FX2);
end
waf(M1+1,:) = S1(1:N).*conj(S2(1:N));
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';


% Normalization
SP1 = fft(hilbert(s1)); 
SP2 = fft(hilbert(s2)); 
indmin = 1+round(fmin*(tcol-2));
indmax = 1+round(fmax*(tcol-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(X)),t,'tfrbert',N,f);
end;