function tfrview(tfr,sig,t,method,param,p1,p2,p3,p4,p5);%TFRVIEW Visualization of time-frequency representations.%        TFRVIEW(TFR,SIG,T,METHOD,PARAM,P1,P2,P3,P4,P5) %        allows to visualize a time-frequency representation.%        TFRVIEW is called through TFRQVIEW from any TFR* function.%%        TFR    : time-frequency representation.%        SIG    : signal in the time-domain.  %        T      : time instants.%        METHOD : chosen representation (name of the corresponding M-file)%        PARAM  : visualization parameter vector :%         PARAM  = [DISPLAY LINLOG THRESHOLD LEVNUMB NF2 LAYOUT FS ISGRID] where%         DISPLAY=1..5 for contour, imagesc, pcolor, surf or mesh%         LINLOG =0/1 for linearly/logarithmically spaced levels%                  %         THRESHOLD  is the visualization threshold, in % %         LEVELNUMB  is the number of levels used with contour%         NF2        is the number of frequency bins displayed%         LAYOUT determines the layout of the figure : TFR alone (1),  %                TFR and SIG (2), TFR and spectrum (3), TFR and SIG and %                spectrum (4), add 4 if you want a colorbar%         FS     is the sampling frequency (may be set to 1.0)%         ISGRID depends on the grids' presence : %                isgrid=isgridsig+2*isgridspec+4*isgridtfr%                where isgridsig=1 if a grid is present on the signal%                and =0 if not, and so on%         fmin   smallest normalized frequency%         fmax   highest  normalized frequency%        P1..P5: parameters of the representation. Run the file %                TFRPARAM(METHOD) to know the meaning of P1..P5.  		 %%	TFRVIEW is called through TFRQVIEW by any file TFR*.%%			Use TFRQVIEW preferably.%			------------------------%%	See also TFRQVIEW, TFRPARAM.%	F. Auger, July 1994, July 1995 - %	O. Lemoine, October-November 1995, May-June 1996. %       F. Auger, May 1998.%	Copyright (c) CNRS - France 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  USAcomp=computer; % so as to know the running computerMatlabVersion=version; MatlabVersion=str2num(MatlabVersion(1));% I hope that future Matlab versions will be more compatibleif (MatlabVersion==4), TickLabelStr='TickLabels';elseif (MatlabVersion>=5), TickLabelStr='Ticklabel';else error('unsupported matlab version. please send an email.'); end; if ( nargin < 5 ), error ('at least 5 parameters are required');end;[tfrrow,tfrcol] = size(tfr); % the size of tfr[trow,tcol]     = size(t);   % the size of t[Nsig,Ncol]     = size(sig); % the size of sigmaxi=max(max(tfr)); % Extraction of the elements of paramdisplay     = param( 1);      % contour, imagesc, pcolor, surf, or meshlinlog      = param( 2);      % linear or logarithmic scalethreshold   = param( 3);      % visualization thresholdlevelnumb   = param( 4);      % number of levelsNf2         = param( 5);      % number of frequency pointslayout      = param( 6);      % figure layoutfs          = param( 7);      % sampling frequencyisgrid      = param( 8);      % grid(s)fmin        = param( 9);      % smallest displayed frequency fmax        = param(10);      % highest displayed frequencyif (fmin>=fmax), error('fmin should be lower than fmax');elseif (fmax>0.5), error('fmax is a normalized frequency, and should be lower than 0.5');endif fs<1000, unitHz=1;                            % display in  s and  Hzelseif (fs>=1e3 & fs<1e6), fs=fs/1e3; unitHz=2;  % display in ms and kHzelseif (fs>=1e6), fs=fs/1e6; unitHz=3;           % display in us and MHzendlinlogtfr=rem(linlog,2);linlogspec=rem(linlog-linlogtfr,4)/2;sigenveloppe=rem(linlog-linlogtfr-linlogspec*2,8)/4;issig=rem(layout-1,2);isspec=rem(layout-1-issig,4)/2;iscolorbar=rem(layout-1-issig-isspec*2,8)/4;if isempty(sig),                        % I can't make miracles issig=0; isspec=0; layout=issig+isspec*2+1;else   layout=layout-4*iscolorbar;end;isgridsig =rem(isgrid,2);isgridspec=rem(isgrid-isgridsig,4)/2;isgridtfr =rem(isgrid-isgridsig-isgridspec*2,8)/4;% Computation of isaffine and freq (vector of frequency samples) if istfraff(method), freq=eval(['p',num2str(nargin-5)]);   	% last input argument is freqs. isaffine=1; if display==2,	                        % imagesc not allowed  display=3;                            % non linear scales for the axes.  disp('Imagesc does not support non-linear scales for axes. We use pcolor instead'); endelse isaffine=0;                            % weyl-heisenberg group of distributions freq=(0.5*(0:Nf2-1)/Nf2);              % dispaly only the positive frequenciesendfreqr=freq*fs; ts=t/fs;                 % real time and frequency% Update variables mini, levels, LinLogStr according to linlog if ~linlogtfr, if (display==4)|(display==5),          % surf or mesh  mini=min(min(tfr)); else  mini=max(min(min(tfr)),maxi*threshold/100.0); end levels=linspace(mini,maxi,levelnumb+1);  LinLogStr=', lin. scale';else mini=max(min(min(tfr)),maxi*threshold/100.0); levels=logspace(log10(mini),log10(maxi),levelnumb+1); LinLogStr=', log. scale';end;% Initialization of the current figurezoom off; clf; set(gcf,'Resize','On','NextPlot','Add');% Initialization of the axesif iscolorbar, axcb   = axes('Units','normal','Visible','off','Box','On'); set(gcf,'UserData',[get(gcf,'UserData') axcb]);end;if issig, axsig  = axes('Units','normal','Visible','off','Box','On'); if comp(1:2)=='PC', set(axsig ,'fontsize',10); end set(gcf,'UserData',[get(gcf,'UserData') axsig]);end;if isspec, axspec = axes('Units','normal','Visible','off','Box','On'); if comp(1:2)=='PC', set(axspec,'fontsize',10); end; set(gcf,'UserData',[get(gcf,'UserData') axspec]);end;axtfr  = axes('Units','normal','Visible','off','Box','On');if comp(1:2)=='PC', set(axtfr ,'fontsize',10); endset(gcf,'UserData',[get(gcf,'UserData') axtfr]); % Test of analycity and computation of spec if ~isempty(sig),  for k=1:Ncol,   isana=1; alpha=2; Lt=max(t)-min(t)+1;   	   if 2*Nf2>=Lt,    spec(:,k)=abs(fft(sig(min(t):max(t),k),2*Nf2)).^2;    else    % modifications :  F. Auger (fog), 30/11/97    Nb_tranches_fog = floor(Lt/(2*Nf2));    % fprintf('%f \n',Nb_tranches_fog);    spec(:,k)=zeros(2*Nf2,1);    for Num_tranche_fog=0:Nb_tranches_fog-1,     % fprintf('%f \n',Num_tranche_fog);     spec(:,k)=spec(:,k)+abs(fft(sig(min(t)+2*Nf2*Num_tranche_fog+(0:2*Nf2-1),k))).^2;    end;    if (Lt>Nb_tranches_fog*2*Nf2),     spectre_fog=fft(sig(min(t)+tfrrow*Nb_tranches_fog:max(t),k),2*Nf2);     spectre_fog=spectre_fog(:);     spec(:,k)=spec(:,k)+abs(spectre_fog).^2;     end;        % sp=abs(fft(sig(min(t):max(t),k))).^2;     % fr1=(0.5*(0:Lt-1)/Lt)*fs;    % fr2=(0.5*(0:2*Nf2-1)/2/Nf2)*fs;    % spec(:,k)=interp1(fr1,sp,fr2);   end   spec1=sum(spec(1:Nf2,k));   spec2=sum(spec(Nf2+1:2*Nf2,k));   if spec2>0.1*spec1,    isana=0;    if ~isreal(sig(min(t):max(t),k)),     alpha=1;    end   end  end end if layout==1,                          % Time-Frequency Representation only  set(axtfr,'Position',[0.10 0.10 0.80 0.80]); elseif layout==2,			% TFR + Signal  set(axtfr,'Position',[0.10 0.10 0.80 0.55]);  set(axsig,'Position',[0.10 0.73 0.80 0.20]);   axes(axsig);    if sigenveloppe,   plot((min(t):max(t))/fs,real(sig(min(t):max(t),:)),...        (min(t):max(t))/fs, abs(sig(min(t):max(t),:)));  else   plot((min(t):max(t))/fs,real(sig(min(t):max(t),:)));  end;  set(gca,['X' TickLabelStr],[]);  ylabel('Real part');  title('Signal in time');  Min=min(min(real(sig))); Max=max(max(real(sig)));  axis([min(ts) max(ts) Min Max]); elseif layout==3,			% TFR + spectrum  set(axspec,'Position',[0.10 0.10 0.15 0.80]);  set(axtfr ,'Position',[0.35 0.10 0.55 0.80]);  axes(axspec);   if isaffine,    f1=freqr(1); f2=freqr(Nf2); df=f2-f1;    Nf4=round((Nf2-1)*fs/(2*df))+1;   for k=1:Ncol,    spec(1:alpha*Nf4,k)=abs(fft(sig(min(t):max(t),k),alpha*Nf4)).^2;    end   spec=spec((round(f1*2*(Nf4-1)/fs)+1):(round(f1*2*(Nf4-1)/fs)+Nf2),:);   freqs=linspace(f1,f2,Nf2);  else   freqs=freqr;   spec=spec(1:Nf2,:);  end  Maxsp=max(max(spec));  if linlogspec==0,   plot(freqs,spec);   title('Linear scale');    % set(axspec,'ytick',[0 Maxsp*max(eps,threshold)*0.01 Maxsp]);   set(axspec,'YTickMode', 'auto');   set(axspec,'Ylim', [Maxsp*threshold*0.01 Maxsp*1.2]);   set(axspec,'Xlim', [fmin fmax]);  else   plot(freqs,10*log10(spec/Maxsp));    title('Log. scale [dB]');   set(axspec,'YTickMode', 'auto');   set(axspec,'Ylim',[10*log10(threshold*0.01) 0]);   set(axspec,'Xlim', [fmin fmax]);  end  xlabel('Energy spectral density');  Nsp=length(spec);   set(gca, ['X' TickLabelStr],[],'view',[-90 90]); elseif layout==4,			% TFR + signal + spectrum  set(axspec,'Position',[0.10 0.10 0.15 0.55]);  set(axsig ,'Position',[0.35 0.73 0.55 0.20]);  set(axtfr ,'Position',[0.35 0.10 0.55 0.55]);  axes(axsig);   if sigenveloppe,   plot((min(t):max(t))/fs,real(sig(min(t):max(t),:)),...        (min(t):max(t))/fs, abs(sig(min(t):max(t),:)));  else   plot((min(t):max(t))/fs,real(sig(min(t):max(t),:)));  end;  ylabel('Real part');  title('Signal in time');  set(gca,['X' TickLabelStr],[]);  Min=min(min(real(sig))); Max=max(max(real(sig)));  axis([min(ts) max(ts) Min Max]);  axes(axspec);   if isaffine, % IF YOU UNDERSTAND SOMETHING TO THAT, PLEASE EXPLAIN ME (f. auger)%   f1=freqr(1); f2=freqr(Nf2); df=f2-f1; %   Nf4=round((Nf2-1)*fs/(2*df))+1;%   for k=1:Ncol,%    spec(1:alpha*Nf4,k)=abs(fft(sig(min(t):max(t),k),alpha*Nf4)).^2; %   end