function [waf,tau,theta] = ambifuwb(X,fmin,fmax,N,trace);%AMBIFUWB Wide-band ambiguity function.%       [WAF,TAU,THETA]=AMBIFUWB(X,FMIN,FMAX,N,TRACE) calculates%	the asymetric wide-band ambiguity function.%%	X     : signal (in time) to be analyzed (the analytic associated%	   signal is considered), of length Nx.%       FMIN,FMAX : respectively lower and upper frequency bounds of%	   the analyzed signal. When specified, these parameters fix%	   the equivalent frequency bandwidth (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 Mellin points (default : automatically determined).%	TRACE : if nonzero, the progression of the algorithm is shown%					(default : 0).%	AF    : matrix containing the coefficients of the ambiguity%	   function. X-coordinate corresponds to the dual variable of %	   scale parameter ; Y-coordinate corresponds to time delay,%	   dual variable of frequency.%	   When called without output arguments, AMBIFUWB displays%	   the squared modulus of the ambiguity function by means of%	   contour.%       TAU   : X-coordinate corresponding to time delay%       THETA : Y-coordinate corresponding to the log(scale) variable%%%	Example :%        sig=altes(128,0.1,0.45); ambifuwb(sig);%%	See also AMBIFUNB.%	P. Goncalvez - December 1995, O. Lemoine - August 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 %	f.auger@ieee.org if (nargin == 0), error('At least one parameter required');end;[Nx,xcol] = size(X);if (xcol==0)|(xcol>1), error('X must have one column');end;if nargin<=4, trace=0; endif (nargin==2), disp('FMIN will not be taken into account. Determine it with FMAX'); disp('     from the following plot of the spectrum.'); elseif (nargin==3), N=[];end;s = hilbert(real(X)); M = round(Nx/2);t = (1:Nx)-M-1;Tmin = 1;Tmax = Nx;T = Tmax-Tmin;if nargin<=2,				        % fmin,fmax,N unspecified STF = fft(fftshift(s));  sp  = (abs(STF(1:M))).^2; Maxsp=max(sp); f   = linspace(0,0.5,M+1) ; f=f(1:M); plot(f,sp) ; grid; xlabel('Normalized frequency'); title('Analyzed signal energy spectrum'); axis([0 1/2 0 1.2*Maxsp]) ;  indmin=min(find(sp>Maxsp/100)); indmax=max(find(sp>Maxsp/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; 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) ; 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 ~isempty(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));geo_f  = fmin*(exp((k-1).*log(q)));tfmatx = zeros(Nx,N);tfmatx = exp(-2*i*t'*geo_f*pi);S = s.'*tfmatx; S = S(ones(1,Nx),:);Sb = S.*tfmatx ;tau = t;S(:,N+1:2*N) = zeros(Nx,N);  S = S.';Sb(:,N+1:2*N) = zeros(Nx,N); Sb = Sb.';% Mellin transform computation of the analyzed signalp=0:(2*N-1);coef = exp(p/2.*log(q))'; MellinS = fftshift(ifft(S(:,1).*coef)).';MellinS = MellinS(ones(1,Nx),:) ; MellinS = MellinS.';for b=1:Nx, if trace, disprog(b,Nx,10); end MellinSb(:,b) = fftshift(ifft(Sb(:,b).*coef)) ;endk = 1:2*N;beta = (p/N-1)/(2*log(q));Scale = logspace(log10(fmin/fmax),log10(fmax/fmin),N) ;waf = zeros(N,Nx) ;MellinSSb = MellinS.*conj(MellinSb) ;waf = ifft(MellinSSb,N);No2=(N+rem(N,2))/2;waf = [waf(No2+1:N,:) ; waf(1:No2,:)]; % NormalizationSP = fft(hilbert(s)); indmin = 1+round(fmin*(Nx-2));indmax = 1+round(fmax*(Nx-2));SPana = SP(indmin:indmax);waf=waf*norm(SPana)^2/waf(No2,M)/N;theta = log(Scale) ;if (nargout==0), contour(tau,theta,abs(waf.^2));  grid xlabel('Delay'); ylabel('Log(scale)'); shading interp title('Wide-band ambiguity function');end;