resgram.m

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

function []=resgram(f,varargin)
%RESGRAM  Reassigned spectrogram plot
%   Usage: resgram(f,op1,op2, ... );
%          resgram(f,sr,op1,op2, ... );
%
%   RESGRAM(f) plots a reassigned spectrogram of f.
%
%   RESGRAM(f,sr) does the same for a signal with sampling rate sr (sampled
%   with sr samples per second);
%
%   Because reassigned spectrograms can be very sparse, consider using the
%   'clim' option (see below).
%
%   Additional arguments can be supplied like this:
%   RESGRAM(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.
%
%   'aa'      - Anti-alias the plot.
%
%   'noaa'    - no anti-aliasing, this is the default.
%
%   '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.
%
%   'lin'     - Show the energy of the coefficients on a linear scale.
%
%   '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
doaa=0;      % Anti-alias.

% 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 'noaa'     
      case 'aa'
        doaa=1;
      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); 

if doeven
  error('Not handled');
else
end;

[itime,ifreq,c]=insttf_dgt(f,{'gauss',tfr_mul},a,M);          
if doaa
  coef=reassign(abs(c).^2,itime,ifreq,a,'aa');
else
  coef=reassign(abs(c).^2,itime,ifreq,a);
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);

  % Remove everything below -100 db of the maximum value.
  %level=30;
  %maxcoef=max(coef(:));
  %coef(coef<maxcoef-level)=maxcoef-level;
  
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;