sgram.m

linear time－frequency toolbox

function []=sgram(f,varargin)
%SGRAM  Spectrogram.
%   Usage: sgram(f,op1,op2, ... );
%
%   SGRAM(f) plots a spectrogram of f using a DGT.
%
%   Additional arguments can be supplied like this:
%   SGRAM(f,'nf','log'). The arguments are character strings,
%   and the available options (so far) are:
%
%-  'nf'      - Display negative frequencies, with the zero-frequency
%               centered in the middle. For real signals, this will just
%               mirror the upper half plane. This is standard for complex
%               signals.
%
%-  'tc'      - Time centering. Move the beginning of the signal to the
%               middle of the plot. This is useful for visualizing the
%               window functions of the toolbox.
%
%-  'db'      - Apply 20*log10 to the coefficients. This makes it possible to
%               see very weak phenomena, but it might show to much noise. A
%               logarithmic scale is more adapted to perception of sound.
%               This is the default of Matlab's SPECGRAM.
%
%-  'vgg'     - Use the signal itself as window. In LaTeX, this is the 
%               symbol V_{g}g This is a quadratic time-frequency
%               distribution related to the Wigner-Ville distribution. 
%
%-  'contour' - Do a contour plot instead of an image.
%          
%-  'surf'    - Do a surf plot instead of an image.
%-
%   In Octave, the default colormap is greyscale. Change it to colormap(jet)
%   for something prettier.

if nargin<1
  error('Too few input arguments.');
end;

% Standard values controlled by optional arguments
dolog=0;
if isreal(f)
  donf=0;
else
  donf=1;
end;

docontour=0;
dosurf=0;
tfr_mul=1;
dotc=0;
dovgg=0;


% Parse optional arguments
if ~isempty(varargin)

  for ii=1:length(varargin)

    optarg=varargin{ii};

    if ischar(optarg)
      switch lower(optarg)
	case 'db'
	  dolog=1;
	case 'nf'
	  donf=1;
	case 'tc'
	  dotc=1;
	  f=fftshift(f);
	case 'vgg'
	  dovgg=1;
	case 'contour'
	  docontour=1;
	case 'surf'
	  dosurf=1;	  
	otherwise
	  error([optarg,' : Unknown optional argument 1']);
      end;
    end;

    if iscell(optarg)
      if isempty(optarg) || ~ischar(optarg{1})
	error(['First element of optinal argument cell array must be a character string.']);
      end;

      switch lower(optarg{1})
	case 'tfr'
	  tfr_mul=optarg{2};
	otherwise
	  error([optarg{1},' : Unknown optional argument 2']);
      end;

    end;
  end;
  
end;

% Approximate resolution along time-axis.
xres=800;

% Ratio of frequency resolution versus time resolution
ratio=3/4;

yres=ceil(xres*ratio);

Ls=length(f);

% Determine time-shift parameter, at least 1.
a=ceil(Ls/xres);

% Number of columns to display
Ndisp=ceil(Ls/a); 

% Determine frequency-shift parameter, at least 1.
b=ceil(Ls/yres);

% Determine transform length.
lcmab=lcm(a,b);

L=ceil(Ls/lcmab)*lcmab;
M=L/b;
N=L/a;

if dovgg
  g=f;
else
  g=pgauss(L,a/b*tfr_mul);
end;

coef=abs(dgt(f,g,a,M)).^2;

if donf
  % Calculate negative frequencies, use DGT
  % Display zero frequency in middle of plot.
 
  % Move zero frequency to the center.
  coef=fftshift(coef,1);
  yr=-L/2:b:L/2;

else
  % Dont calculate negative frequencies, try to use DGT of twice the size.
  coef=coef(1:M/2,:);
  yr=0:b:L/2;

end;

if dotc
  xr=-floor(Ndisp/2)*a:a:floor((Ndisp-1)/2)*a;
else
  xr=0:a:Ndisp*a-1;
end;

% Cut away zero-extension.
coef=coef(:,1:Ndisp);

% Apply transformation to coefficients.
if dolog
  % We have already squared the coefficients, so only multiply
  % by 10. We add eps to avoid taking the log of zero.
  coef=10*log10(coef+eps);
end;

if isoctave
  % Do Octave plotting.
  imagesc(flipud(coef));
else
  % Do Matlab plotting.

  if docontour
    % XXX Needs xr and yr
    ih=contour(coef);
  else
    if dosurf
      % XXX Needs xr and yr
      ih=surf(coef);
    else      
      ih=imagesc(xr,yr,coef);
    end;
  end;
  axis('xy');
  xlabel('Time')
  ylabel('Frequency')

end;