sgram.m

Matlab时频分析工具箱,希望能对大家有所帮助啊

function []=sgram(f,varargin)
%SGRAM  Spectrogram.
%   Usage: sgram(f,op1,op2, ... );
%          sgram(f,sr,op1,op2, ... );
%
%   SGRAM(f) plots a spectrogram of f using a DGT.
%
%   SGRAM(f,sr) does the same for a signal with sampling rate sr (sampled
%   with sr samples per second);
%
%   Additional arguments can be supplied like this:
%   SGRAM(f,'nf','tfr',2,'log'). The arguments must be character strings
%   possibly followed by an argument:
%
%   'tfr',v   - Set the ratio of frequency resolution to time resolution.
%               A value v=1 is the default. Setting v>1 will give better
%               frequency resolution at the expense of a worse time
%               resolution. A value of 0<v<1 will do the opposite.
%
%   'thr',r   - Keep only the largest fraction r of the coefficients, and
%               set the rest to zero.
%
%   'p'       - Use a periodic bounday condition. This is the default.
%
%   'e'       - Use an even boundary condition. This produces the fewest
%               visible artifacts at the boundary. Note: Currently this
%               doubles the computational time.
%
%   '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.
%
%   'lin'     - Show the energy of the coefficients on a linear scale.
%
%   '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. 
%
%   'image'   - Use 'imagesc' to display the spectrogram. This is the
%               default.
%
%   'clim',[clow,chigh] - Use a colormap ranging from clow to chigh. These
%               values are passed to IMAGESC. See the help on IMAGESC.
%
%   'range',range - Use a colormap in the interval [chigh-range,chigh], where
%               chigh is the highest value in the plot.
%
%   'fmax',y  - Display y as the highest frequency.
%
%   'contour' - Do a contour plot to display the spectrogram.
%          
%   'surf'    - Do a surf plot to display the spectrogram.
%
% 
%   In Octave, the default colormap is greyscale. Change it to colormap(jet)
%   for something prettier.

% 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 3 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, see <http://www.gnu.org/licenses/>.

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

if sum(size(f)>1)>1
  error('Input must be a vector.');
end;

global TF_CONF;
% Approximate resolution along time-axis.
xres=TF_CONF.xres;

% Ratio of frequency resolution versus time resolution
displayratio=TF_CONF.displayratio;

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

dolog=1;     % Use a Db scale.
docontour=0; % Do a contour plot.
dosurf=0;    % Do a surf plot.
dotc=0;      % Display the beginning of the signal in the middle.
dovgg=0;     % Use the signal as window.
dosr=0;      % Display sampling rate
domask=0;    % Remove small coefficients.
doeven=0;    % Even extension of signal to reduce boundary effect.
doclim=0;    % Limit the colorscale, see IMAGESC for help.
dorange=0;   % Limit the colorscale.
dofmax=0;    % Limit the frequency axis

% Set the default time/frequency ratio.
tfr_mul=1;   

% Resampling rate: Used when fmax is issued.
resamp=1;

% Parse optional arguments
if ~isempty(varargin)

  startii=1;

  if isnumeric(varargin{1})
    dosr=1;
    sr=varargin{1};
    startii=2;
  end;
  
  ii=startii-1;
  while ii<length(varargin)   

    ii=ii+1;

    optarg=varargin{ii};
    
    if ~ischar(optarg)
       error('Argument %i must be a character string.',ii+1);
    end;
    
    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;	  
      case 'lin'
	dolog=0;
      case 'thr'
	domask=1;
      case 'image'
      case 'p'
      case 'e'
	doeven=1;
      case 'tfr'
	tfr_mul=varargin{ii+1};
	ii=ii+1;
      case 'thr'
	domask=1;
	mask_val=varargin{ii+1};
	ii=ii+1;
      case 'clim'
	doclim=1;
	clim=varargin{ii+1};
	ii=ii+1;
      case 'range'
	dorange=1;
	range=varargin{ii+1};
	ii=ii+1;
      case 'fmax'
	dofmax=1;
	fmax=varargin{ii+1};
	ii=ii+1;
      otherwise
	error([optarg,' : Unknown optional argument 1']);
    end;
 
 end;
  
end;

yres=ceil(xres*displayratio);

% Downsample
if dofmax
  if dosr
    resamp=fmax*2/sr;
  else
    resamp=fmax*2/length(f);
  end;

  f=fftresample(f,round(length(f)*resamp));
end;

Ls=length(f);

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

% Determine number of channels.
M=min(yres,Ls);

L=iirpar(Ls,a,M);
N=L/a;
b=L/M;

% 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;

if doeven
  if dovgg
    error('XXX');
  else
    if mod(a,2)==0
      g=pgauss(2*L,.5*tfr_mul,.5);
    else
      g=pgauss(2*L,.5*tfr_mul,0);
    end;
  end;
else
  if dovgg
    g=f;
  else
    g=pgauss(L,tfr_mul);
  end;
end;

if doeven
  f=[f;flipud(f)];
  eL=length(f);
  f=f.*(exp(-2*pi*i*(eL-1)).');
  g=circshift(g,floor(a/2));
  coef=abs(dgt(f,g,a,M)).^2;
  % Phase correct is not included because we only want the magnitude.
  %coef=coef(:,1:size(coef,2)/2);
else
  coef=abs(dgt(f,g,a,M)).^2;
end;

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

  if dofmax
    yr=-fmax:fmax/M:fmax;
  else
    if dosr
      yr=-sr/2:sr/M:sr/2;
    else
      yr=-L/2:b:L/2;
    end;
  end;
else
  % Dont calculate negative frequencies, try to use DGT of twice the size.
  coef=coef(1:ceil((M+1)/2),:);
  if dofmax
    yr=0:fmax/M:fmax;
  else
    if dosr
      yr=0:sr/M:sr/2;
    else
      yr=0:b:L/2;
    end;
  end;
end;

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

if dosr
  % Scale x-axis by sampling rate.
  xr=xr/sr;  
end;

% Scale x-axis by resampling rate
xr=xr/resamp;

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

if domask
  % keep only the largest coefficients.
  coef=largestr(coef,mask_val);
end

% 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+realmin);
end;

% 'range' parameter is handled by converting it to clim.
if dorange
  maxclim=max(coef(:));
  clim=[maxclim-range,maxclim];
  doclim=1;
end;

if isoctave
  % Do Octave plotting.
  if doclim
    imagesc(flipud(coef),clim);
  else
    imagesc(flipud(coef));
  end;
else
  % Do Matlab plotting.

  if docontour
    % XXX Needs xr and yr
    contour(coef);
  else
    if dosurf
      % XXX Needs xr and yr
      surf(coef);
    else      
      if doclim
	imagesc(xr,yr,coef,clim);
      else
	imagesc(xr,yr,coef);
      end;
    end;
  end;
  axis('xy');
  if dosr
    xlabel('Time (s)')
    ylabel('Frequency (Hz)')
  else
    xlabel('Time')
    ylabel('Frequency')
  end;

end;