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 %	f.auger@ieee.org global ratio_fif (nargin == 0), error('At least one parameter required');end;[xrow,xcol] = size(X);if nargin<=6, trace=0; endif (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);ends1 = 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 fmin==[], fmin=fmindflt; end if fmax==[], fmax=fmaxdflt; endendif 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');endB = 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);endif N~=[], if (N<Nmin),  dispstr=['Warning : the number of analyzed voices (N) should be > ',num2str(Nmin)];  disp(dispstr); endelse N=Ndflt; endfmin_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 spectrumk = 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 tfmatxS1(N+1:2*N) = zeros(1,N);S2(N+1:2*N) = zeros(1,N);% Mellin transform computation of the analyzed signalp = 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 signalwaf = 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));  endendwaf  = [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 interpolationtfr  = 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); tfr(n,:) = interp1(x,y,xi,'spline')';end t = time;f = geo_f(1:N)';% NormalizationSP1 = 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_fif (nargout==0), tfrqview(tfr,hilbert(real(X1)),t,'TFRUNTER',N,f);end;