tfrview.m

最新时频分析处理软件

function tfrview(tfr,sig,t,method,param,map,p1,p2,p3,p4,p5);
%TFRVIEW Visualization of time-frequency representations.
%	TFRVIEW(TFR,SIG,T,METHOD,PARAM,MAP,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 
%				ACCESS STATE 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/remove the colorbar (5)
%	  ACCESS depends on the way you access to tfrview : from the 
%		command line (0) ; from tfrqview, except after a
%		change in the sampling frequency or in the layout (1) ; 
%		from tfrqview, after a change in the layout (2) ; 
%		from tfrqview, after a change in the sampling frequency (3)
%	  STATE depends on the signal/colorbar presence : 
%		no signal, no colorbar (0) ; signal, no colorbar (1) ;
%		no signal, colorbar (2) ; signal and colorbar (3)
%	  FS is the sampling frequency (may be set to 1.0)
%	  ISGRID depends on the grids' presence : 
%			isgrid=isgridsig+2*isgridspe+4*isgridtfr
%		where isgridsig=1 if a grid is present on the signal
%		and =0 if not, and so on
%	MAP    :  selected colormap. 
%	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. 
%	Copyright (c) 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 
%	[e-mail] 

comp=computer;
 
if ( nargin < 6 ),
 error ('at least 6 parameters are required');
end;
[tfrrow,tfrcol] = size(tfr);
[trow,tcol] = size(t);
[Nsig,Ncol]=size(sig);
maxi=max(max(tfr)); 

% Extraction of the elements of param
display   = param(1); 
linlog    = param(2); 
threshold = param(3); 
levelnumb = param(4);
Nf2       = param(5);
layout	  = param(6);
access    = param(7);
state     = param(8);
fs        = param(9);
isgrid    = param(10);

if fs<1000,
 unitHz=1;
elseif (fs>=1e3 & fs<1e6),
 fs=fs/1e3;
 unitHz=2;
elseif (fs>=1e6),
 fs=fs/1e6;
 unitHz=3;
end

linlogtfr=rem(linlog,2);
linlogspec=(linlog-linlogtfr)/2;
issig=rem(state,2);
iscolorbar=(state-issig)/2;
isgridsig=rem(isgrid,2);
tempo=(isgrid-isgridsig)/2;
isgridspe=rem(tempo,2);
isgridtfr=(tempo-isgridspe)/2;


% Computation of isaffine and freq (vector of frequency samples) 
if strcmp(method,'TFRASPW' ) | strcmp(method,'TFRSCALO') | ...
   strcmp(method,'TFRDFLA' ) | strcmp(method,'TFRSPAW' ) | ...
   strcmp(method,'TFRUNTER') | strcmp(method,'TFRBERT' ),
 freq=eval(['p',num2str(nargin-6)]);   	% last input argument is freqs.
 isaffine=1;
 if display==2,				% imagesc do not allow
  display=3;				% non linear scales for the axes.
  disp('Imagesc does not support non-linear scales for axes. We use pcolor instead');
 end
else
 isaffine=0;
 freq=(0.5*(0:Nf2-1)/Nf2);
end
freqr=freq*fs;
ts=t/fs;


% Update variables mini, levels, Strlinlog according to linlog 
if ~linlogtfr,
 if (display==4|display==5),
  mini=min(min(tfr));
 else
  mini=max(min(min(tfr)),maxi*threshold/100.0);
 end
 levels=linspace(mini,maxi,levelnumb+1); 
 Strlinlog=', lin. scale';
else
 mini=max(min(min(tfr)),maxi*threshold/100.0);
 levels=logspace(log10(mini),log10(maxi),levelnumb+1);
 Strlinlog=', log. scale';
end;


% Initialization of the axes
if (access==1|access==2|access==3),			% From the main menu of tfrqview
 axesh  = get(gcf,'UserData');
 axcb   = axesh(1);
 axsig  = axesh(2);
 axspec = axesh(3);
 axtfr  = axesh(4);
elseif (access==0),		% From the save option of tfrqview
 axcb   = axes('Units','normal','Visible','off','Box','On');
 axsig  = axes('Units','normal','Visible','off','Box','On');
 axspec = axes('Units','normal','Visible','off','Box','On');
 axtfr  = axes('Units','normal','Visible','off','Box','On');
end

if (access==0|access==2|access==3),

 % Test of analycity and computation of spec
 if 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
    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>spec1/10,
    isana=0;
    if ~isreal(sig(min(t):max(t),k)),
     alpha=1;
    end
   end
  end
 end

 if layout==1,			% Time-Frequency Representation

  set(axtfr,'Position',[0.13 0.11 0.775 0.815]);
  issig=0; isspec=0;

 elseif layout==2,			% TFR + Signal

  set(axtfr,'Position',[.1 .1 .83 .55]);
  set(axsig,'Position',[.1 .72 .83 .2]); 
  axes(axsig);  
  plot(ts,real(sig(t,:)));
  set(gca,'xticklabels',[]);
  if comp(1:2)=='PC', set(gca,'fontsize',10); end
  ylabel('Real part');
  title('Signal in time');
  Min=min(min(real(sig))); Max=max(max(real(sig)));
  axis([min(ts) max(ts) Min Max]);
  issig=1; isspec=0;

 elseif layout==3,			% TFR + spectrum

  axes(axspec); 
  set(gca,'Position',[.1 .1 .18 .8]);
  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,
   plot(freqs,spec);
   if comp(1:2)=='PC', set(gca,'fontsize',10); end
   title('Linear scale'); 
   set(gca,'ytick',[0 round(Maxsp/2) round(Maxsp)]);
  else
   semilogy(freqs,spec);
   if comp(1:2)=='PC', set(gca,'fontsize',10); end
   title('Log. scale [dB]');
   set(gca,'ytick',[fix(10^(log10(Maxsp)/2));fix(Maxsp)]);
   str1=['   ']; str2=['   '];
   st1=num2str(fix(10*log10(Maxsp)/2)); str1(1:length(st1))=st1;
   st2=num2str(fix(10*log10(Maxsp)));   str2(1:length(st2))=st2; 
   set(axspec,'yticklabels',[str1;str2]);
  end
  xlabel('Energy spectral density');
  Nsp=length(spec); 
  set(gca,'Xlim',[freqs(1) freqs(Nsp)]);
  set(gca,'Ylim',[0 Maxsp*1.2]);
  set(gca,'xticklabels',[],'view',[-90 90]);

  set(axtfr,'Position',[.36 .1 .57 .8]);
  issig=0; isspec=1; 

 elseif layout==4,			% TFR + signal + spectrum
 
  set(axsig,'Position',[.33 .73 .6 .2]);
  axes(axsig); 
  plot(ts,real(sig(t,:)));
  if comp(1:2)=='PC', set(gca,'fontsize',10); end
  ylabel('Real part');
  title('Signal in time');
  set(gca,'xticklabels',[]);
  Min=min(min(real(sig))); Max=max(max(real(sig)));
  axis([min(ts) max(ts) Min Max]);

  set(axspec,'Position',[.1 .1 .15 .55]);
  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,
   plot(freqs,spec);
   if comp(1:2)=='PC', set(gca,'fontsize',10); end
   title('Linear scale');  
   set(axspec,'ytick',[0 round(Maxsp/2) round(Maxsp)]);
  else
   semilogy(freqs,spec); 
   if comp(1:2)=='PC', set(gca,'fontsize',10); end
   title('Log. scale [dB]');
   set(gca,'ytick',[fix(10^(log10(Maxsp)/2));fix(Maxsp)]);
   str1=['   ']; str2=['   '];
   st1=num2str(fix(10*log10(Maxsp)/2)); str1(1:length(st1))=st1;
   st2=num2str(fix(10*log10(Maxsp)));   str2(1:length(st2))=st2; 
   set(axspec,'yticklabels',[str1;str2]);
  end
  xlabel('Energy spectral density');
  Nsp=length(spec); 
  set(gca,'Xlim',[freqs(1) freqs(Nsp)]);
  set(gca,'Ylim',[0 Maxsp*1.2]);
  set(gca,'xticklabels',[],'view',[-90 90]);

  set(axtfr,'Position',[.33 .1 .6 .55]);
  issig=1; isspec=1;

 elseif layout==5,			% Colorbar

  if ~iscolorbar,                     % there is no color bar already 
   pos=get(axtfr,'Position');
   set(axtfr,'Position',[pos(1) pos(2) pos(3)-0.03 pos(4)]);
   set(axcb,'Position',[pos(1)+pos(3)-0.01 pos(2) 0.01 pos(4)]);
   axes(axcb); 
   Ncolors=length(colormap);