phaseplot.m

来自「linear time－frequency toolbox」· M 代码 · 共 208 行

function []=phaseplot(f,varargin)
%PHASEPLOT  Phase plot
%   Usage: phaseplot(f,op1,op2, ... );
%
%   PHASEPLOT(f) plots the phase of f using a DGT.
%
%   Additional arguments can be supplied like this:
%   PHASEPLOT(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 usefull for visualizing the
%               window functions of the toolbox.
%
%-  'vgg'     - Use the signal itself as window. In LaTeX, this is the 
%               symbol V_{g}g
%
%-  'thr'     - Use the amplitude values to threshold the phase values.
%               For small amplitude values the phase values are. This is
%               because an error, that is small in amplitude can still
%               lead to arbitrary phase values in the vicinity of zero.
%               Setting this flag will set the phase is to a constant value
%               (0) for all pairs (m,n), where the amplitude is below a 
%               certain relative value, set to 0.001 of the maximum amplitude.
%
%-  'phl'     - Different to dgt.m the Gabor transform can be defined 
%               by shifting the order of translation and modulation. This
%               is sometimes used in algorithms to be able to directly apply
%               the FFT without circular shifts, as for every time point the FFT
%               can start from zero. Obviously this choice has a big influence
%               on the development of the phase of the coefficients over time. 
%
%-
%   In Octave, the default colormap is greyscale. Change it to colormap(jet)
%   for something prettier.
%
%   Obviously in the spectrogram the information about the complex angle of 
%   the Gabor coefficients, i.e. the phase, gets lost. To plot the phase 
%   some issues have to be considered:
% 
%   - phase unwrapping: The angle function is non-continuous on a half-axis
%   in the complex plane. So there are continuous complex function, whose 
%   angle is not continuous, e.g. a periodic rotation around zero. A 
%   work-around is to add integer multiples of 2 pi to the result, the so 
%   called 'unwrapping'. For plots a continuous effect can just be reached 
%   by choosing a fitting color scheme. 
% 
%   SEE ALSO: PHASELOCK
%
%   EXAMPLES: EXAMP_PHASEPLOT
% 
%   REFERENCES:
%     R. Carmona, W. Hwang, and B. Torrésani. Practical Time-Frequency Analysis:
%     continuous wavelet and Gabor transforms, with an implementation in S,
%     volume 9 of Wavelet Analysis and its Applications. Academic Press, San
%     Diego, 1998.

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

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

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


% Parse optional arguments
if ~isempty(varargin)

  for ii=1:length(varargin)

    optarg=varargin{ii};

    if ischar(optarg)
      switch lower(optarg)
	case 'nf'
	  donf=1;
	case 'tc'
	  dotc=1;
	  f=fftshift(f);
	case 'vgg'
	  dovgg=1;
    case 'thr'
      domask=1;
    case 'phl'
      dophaslock=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;

comp_coef=dgt(f,g,a,M);

if dophaslock
  % use phase locking
  comp_coef = phaselock(comp_coef,a);
end
coef = angle(comp_coef);

if domask
  % use threshold based on amplitude
  ampl_coef = comp_coef.^2;
  mma = max(max(ampl_coef));
  indx = find(ampl_coef <= 1e-3*mma);
  coef(indx)=0;
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);
  yr=-L/2:b:L/2;
else
  % Dont calculate negative frequencies, try to use DGT of twice the size.

  coef=coef(floor(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);

if isoctave
  % Do Octave plotting.
  imagesc(flipud(coef));
else
  % Do Matlab plotting.
  ih=imagesc(xr,yr,coef);
  axis('xy');
  xlabel('Time')
  ylabel('Frequency')

end;