function [tfr,t,f,wt]=tfrscalo(X,time,wave,fmin,fmax,N,trace);%TFRSCALO Scalogram, for Morlet or Mexican hat wavelet.%	[TFR,T,F,WT]=TFRSCALO(X,T,WAVE,FMIN,FMAX,N,TRACE) computes %	the scalogram (squared magnitude of a continuous wavelet%	transform). %%	X : signal (in time) to be analyzed (Nx=length(X)). Its%	    analytic version is used (z=hilbert(real(X))).  %	T : time instant(s) on which the TFR is evaluated %	     					(default : 1:Nx).%	WAVE : half length of the Morlet analyzing wavelet at coarsest % 	    scale. If WAVE = 0, the Mexican hat is used. WAVE can also be%           a vector containing the time samples of any bandpass%           function, at any scale.        	(default : sqrt(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.%	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, TFRSCALO%	    runs TFRQVIEW.%	F : vector of normalized frequencies (geometrically sampled %	    from FMIN to FMAX).%	WT : Complex matrix containing the corresponding wavelet%	    transform. The scalogram TFR is the square modulus of WT.%%	Example :    %	 sig=altes(64,0.1,0.45); tfrscalo(sig);  %%	See also all the time-frequency representations listed in%	the file CONTENTS (TFR*)%	P. Goncalves, October 1995 - O. Lemoine, June 1996. %	Copyright (c) 1995 Rice University - CNRS 1996.%%  This program is free software; you can redistribute it and/or modify%  it under the terms of the GNU General Public License as published by%  the Free Software Foundation; either version 2 of the License, or%  (at your option) any later version.%%  This program is distributed in the hope that it will be useful,%  but WITHOUT ANY WARRANTY; without even the implied warranty of%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the%  GNU General Public License for more details.%%  You should have received a copy of the GNU General Public License%  along with this program; if not, write to the Free Software%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USAif (nargin == 0), error('At least one parameter required');end;[xrow,xcol] = size(X);if nargin<=6, trace=0; endif (nargin == 1), time=1:xrow; wave=sqrt(xrow);elseif (nargin == 2), wave=sqrt(xrow);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=xrow;end;[trow,tcol] = size(time);if (xcol==0)|(xcol>2), error('X must have one or two columns');elseif (trow~=1), error('TIME must only have one row'); elseif wave<0, error('WAVE must be positive');end; s = (real(X) - mean(real(X)))';  z = hilbert(s) ;if trace, disp('Scalogram distribution'); end;if nargin<=4		        % fmin,fmax,N unspecified STF = fft(fftshift(z(min(time):max(time)))); Nstf=length(STF); sp = (abs(STF(1:round(Nstf/2)))).^2; Maxsp=max(sp); f = linspace(0,0.5,round(Nstf/2)+1) ; f = f(1:round(Nstf/2)); plot(f,sp) ; grid; xlabel('Normalized frequency'); title('Analyzed signal energy spectrum'); 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; end txt=['Number of frequency samples [',num2str(2^nextpow2(xrow)),'] : '];  N=input(txt);  if isempty(N), N=2^nextpow2(xrow); endendfmin_s=num2str(fmin); fmax_s=num2str(fmax); N_s=num2str(N);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');endif trace, disp(['Frequency runs from ',fmin_s,' to ',fmax_s,' with ',N_s,' points']);endf = logspace(log10(fmin),log10(fmax),N);a = logspace(log10(fmax/fmin),log10(1),N); wt =zeros(N,tcol);tfr=zeros(N,tcol);if wave > 0 if trace, disp(['using a Morlet wavelet']) ; end for ptr=1:N,  if trace, disprog(ptr,N,10); end  nha = wave*a(ptr);  tha = -round(nha) : round(nha);  ha  = exp(-(2*log(10)/nha^2)*tha.^2).*exp(i*2*pi*f(ptr)*tha);   detail = conv(z,ha)./sqrt(a(ptr));  detail = detail(round(nha)+1:length(detail)-round(nha)) ;  wt(ptr,:)  = detail(time) ;  tfr(ptr,:) = detail(time).*conj(detail(time)) ; endelseif wave == 0 if trace, disp(['using a Mexican hat wavelet']) ; end for ptr = 1:N  if trace, disprog(ptr,N,10); end  ha  = mexhat(f(ptr)) ;  nha = (length(ha)-1)/2 ;  detail = conv(z,ha)./sqrt(a(ptr));  detail = detail(round(nha)+1:length(detail)-round(nha)) ;  wt(ptr,:)  = detail(time);  tfr(ptr,:) = detail(time).*conj(detail(time)) ; end  elseif length(wave) > 1 [rwav,cwav]=size(wave); if cwav>rwav, wave=wave.'; end wavef = fft(wave) ; nwave = length(wave) ; f0 = find(abs(wavef(1:nwave/2)) == max(abs(wavef(1:nwave/2)))); f0 = mean((f0-1).*(1/nwave)); if trace, disp(['mother wavelet centered at f0 = ',num2str(f0)]); end a = logspace(log10(f0/fmin),log10(f0/fmax),N); B = 0.99; R = B/((1.001)/2);  nscale = max(128,round((B*nwave*(1+2/R)*log((1+R/2)/(1-R/2)))/2)); if trace, disp('Scale computation :'); end wts = scale(wave,a,fmin,fmax,nscale,trace); for ptr = 1:N,   clear detail  if trace, disprog(ptr,N,10); end  ha = wts(ptr,:);  nha = length(ha)/2;  detail = conv(z,ha)./sqrt(a(ptr));  detail = detail(fix(nha):length(detail)-round(nha));  wt(ptr,:) = detail(time);  tfr(ptr,:) = detail(time).*conj(detail(time)); endendt = time;f = f';% NormalizationSP = fft(z); indmin = 1+round(fmin*(xrow-2));indmax = 1+round(fmax*(xrow-2));SPana=SP(indmin:indmax);tfr=tfr*norm(SPana)^2/integ2d(tfr,t,f)/N;if (nargout==0), tfrqview(tfr,hilbert(real(X)),t,'tfrscalo',wave,N,f);end;