function [tfr,t,f]=tfrspaw(X,time,K,nh0,ng0,fmin,fmax,N,trace);%TFRSPAW Smoothed Pseudo Affine Wigner time-frequency distributions.%	[TFR,T,F]=TFRSPAW(X,T,K,NH0,NG0,FMIN,FMAX,N,TRACE)%	generates the auto- or cross- Smoothed Pseudo Affine Wigner%	distributions.  %%	X : signal (in time) to be analyzed. If X=[X1 X2], TFRSPAW %	   computes the cross-Smoothed Pseudo Affine Wigner distribution.%						(Nx=length(X)).%	T : time instant(s) on which the TFR is evaluated  %	   					(default : 1:Nx).%	K : label of the K-Bertrand distribution. The distribution with%	   parameterization function %	   lambdak(u,K) = (K (exp(-u)-1)/(exp(-Ku)-1))^(1/(K-1)) %	   is computed				(default : 0).%	     K=-1 : Smoothed pseudo (active) Unterberger distribution %	     K=0  : Smoothed pseudo Bertrand distribution%	     K=1/2: Smoothed pseudo D-Flandrin distribution%	     K=2  : Affine smoothed pseudo Wigner-Ville distribution.%	NH0 : half length of the analyzing wavelet at coarsest scale.  %	   A Morlet wavelet is used. NH0 controles the frequency %	   smoothing of the smoothed pseudo Affine Wigner distribution.%						(default : sqrt(Nx)).%	NG0 : half length of the time smoothing window. %	   NG0 = 0 corresponds to the Pseudo Affine Wigner distribution.  %						(default : 0).%	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, TFRSPAW%	   runs TFRQVIEW.%	F : vector of normalized frequencies (geometrically sampled %	   from FMIN to FMAX).%%	Example :    %	 sig=altes(64,0.1,0.45); tfrspaw(sig);%	 %	See also all the time-frequency representations listed in%	the file CONTENTS (TFR*)%	P. Goncalves, October 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 if (nargin == 0), error('At least one parameter required');end;[xrow,xcol] = size(X);if (nargin<=8), trace=0; endif (nargin == 1), time=1:xrow; K=0; nh0=sqrt(xrow); ng0=0;elseif (nargin == 2), K=0; nh0=sqrt(xrow); ng0=0;elseif (nargin == 3), nh0=sqrt(xrow); ng0=0;elseif (nargin == 4), ng0=0;elseif (nargin == 6), disp('FMIN will not be taken into account. Determine it with FMAX'); disp('     from the following plot of the spectrum.'); elseif (nargin == 7), 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,  if (K==-1),  disp('Smoothed pseudo (active) Unterberger distribution'); elseif K==0,   disp('Smoothed pseudo Bertrand distribution'); elseif K==1/2,   disp('Smoothed pseudo D-Flandrin distribution'); elseif K==2,   disp('Affine smoothed pseudo Wigner-Ville distribution'); else  disp('Smoothed Pseudo Affine Wigner distribution'); end;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;if nargin<=6,				        % 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; 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) ; Qte  = fmax/fmin ;    umax = log(Qte); Teq  = nh0/(fmax*umax);  if Teq<2*nh0, M0 = (2*nh0^2)/Teq-nh0+1;else M0 = 0;end;MU = round(nh0+M0);T  = 2*MU-1;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<=6, 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']);endk = 1:N;q = (fmax/fmin)^(1/(N-1));a = exp((k-1).*log(q));         % a is an increasing scale vector.geo_f = fmin*a;                 % geo_f is a geometrical increasing                                % frequency vector.% Wavelet decomposition computationmatxte1 = zeros(N,tcol);matxte2 = zeros(N,tcol);[p1,p2,p3,wt1] = tfrscalo(s1,time,nh0,fmin,fmax,N) ;[p1,p2,p3,wt2] = tfrscalo(s2,time,nh0,fmin,fmax,N) ;for ptr = 1:N,  matxte1(ptr,:) = wt1(ptr,:).*sqrt(a(N-ptr+1)) ;  matxte2(ptr,:) = wt2(ptr,:).*sqrt(a(N-ptr+1)) ; end ;umin = -umax;u=linspace(umin,umax,2*MU+1);du = u(2)-u(1);u=u(1:2*MU);u(MU+1) = 0;p = 0:(2*N-1);beta = (p/N-1)./(2*log(q));l1=zeros(2*MU,2*N);l2=zeros(2*MU,2*N);for m = 1:2*MU, l1(m,:) = exp(-2*i*pi*beta*log(lambdak( u(m),K))); l2(m,:) = exp(-2*i*pi*beta*log(lambdak(-u(m),K)));end % Determination of the time smoothing window Gif ng0==0, G = ones(2*MU,1);else a_t = 3 ;            % (attenuation of 10^(-a_t) at t = tmax) sigma_t = ng0*fmax/sqrt(2*a_t*log(10)); a_u = 2 * pi^2 * sigma_t^2 * umax^2 / log(10) ; sigma_u = 1/(2 * pi * sigma_t) ; G = exp(-(a_u*log(10)/MU^2)*[-MU:MU-1].^2);  if sigma_u < du  disp('Maximum time smoothing reached. Increase width of wavelet for effectiveness.') ; end G=G';endwaf = zeros(2*MU,N);tfr = zeros(N,tcol);S1  = zeros(1,2*N);S2  = zeros(1,2*N);MX1 = zeros(2*N,2*MU);MX2 = zeros(2*N,2*MU);TX1 = zeros(2*MU,N);TX2 = zeros(2*MU,N);for ti = 1:tcol, if trace, disprog(ti,tcol,10); end S1(1:N) = matxte1(:,ti).'; Mellin1 = fftshift(ifft(S1)); MX1 = (l1.*Mellin1(ones(1,2*MU),:)).'; MX1 = fft(MX1); TX1 = MX1(1:N,:).'; S2(1:N) = matxte2(:,ti).'; Mellin2 = fftshift(ifft(S2));      MX2 = (l2.*Mellin2(ones(1,2*MU),:)).'; MX2 = fft(MX2); TX2 = MX2(1:N,:).'; waf = real(TX1.*conj(TX2)).*G(:,ones(N,1)); tfr(:,ti) = (sum(waf).*geo_f).';	% first row of tfr corresponds to				% the lowest frequency.end;t = time; f = geo_f'; % 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;if (nargout==0), tfrqview(real(tfr),hilbert(real(X)),t,'tfrspaw',K,nh0,ng0,N,f);end;