AUDSPECGRAM - Auditory spectrogram.

Program code:

function varargout=audspecgram(insig,fs,varargin)

%   AUTHOR : Peter Soendergaard.
%   TESTING: NA
%   REFERENCE: NA

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

if ~isnumeric(insig) || ~isvector(insig)
  error('%s: Input must be a vector.',upper(mfilename));
end;

global CASP_CONF;

if isempty(CASP_CONF)
  error(['You need to run CASPSTART. This could be because you have used a ', ...
        'CLEAR ALL statement.']);
end;

% Get the default values.
defaults=CASP_CONF.plotdefaults;

% Switches to turn specific behaviours on and off.
doadapt=1;   % Model adaptation
dorectify=1; % Use half-wave rectification.
dolog=1;     % Use a Db scale.
domask=0;    % Remove small coefficients.
doclim=0;    % Limit the colorscale, see IMAGESC for help.
dorange=0;   % Limit the colorscale.
doxres=0;
doyres=0;
domlp=0;     % Modulation low-pass
dombp=0;     % Modulation band-pass

% Parse default and optional arguments
arguments={defaults{:},varargin{:}};
startii=1;

ii=startii-1;
while ii2*1.5*fhigh

  fsnew=round(fhigh*2*1.2);

  % Determine new signal length
  siglen=round(siglen/fs*fsnew);

  % Do the resampling using an FFT based method, as this is more flexible
  % than the 'resample' method included in Matlab
  insig=fftresample(insig,siglen);

  % Switch to new value
  fs=fsnew;
end;

% Determine the hopsize
% Using a hopsize different from 1 is currently not possible because all
% the subsequent filters fail because of a to low subband sampling rate.
%hopsize=max(1,floor(siglen/xres));

hopsize=1;

% find the center frequencies used in the filterbank
fc = erbspace(flow,fhigh,yres);

% Calculate filter coefficients for the gammatone filter bank.
[gt_b, gt_a]=gammatone(fc, fs);

% Apply the Gammatone filterbank
outsig = filterbank(gt_b,gt_a,insig,hopsize);

% The subband are now (possibly) sampled at a lower frequency than the
% original signal.
fssubband=round(fs/hopsize);

if dorectify && (fssubband>2000)
  outsig = envextract(2*real(outsig),fssubband);
else
  outsig = abs(outsig);
end;

if doadapt
  % non-linear adaptation loops
  outsig = adaptloop(outsig, fssubband,10);
end;

if domlp
  % Calculate filter coefficients for the 50 Hz modulation lowpass filter.
  mlp_a = exp(-mlp/fssubband);
  mlp_b = 1 - mlp_a;
  mlp_a = [1, -mlp_a];

  % Apply the low-pass modulation filter.
  outsig = filter(mlp_b,mlp_a,outsig);
end;

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

% Apply transformation to coefficients.
if dolog
  % This is a safety measure to avoid log of negative numbers if the
  % users chooses an incorrect combination of options (e.g. 'adapt' and 'log').
  outsig(:)=max(outsig(:),realmin);

  outsig=20*log10(outsig);
end;

% 'dynrange' parameter is handled by threshholding the coefficients.
if dorange
  maxclim=max(outsig(:));
  outsig(outsig0
  varargout={outsig,fc};
end;