% gtep: Gammatone filter bank excitation pattern
% Usage: [cf,expat]=gtep(y,rate,nfilt,flo,fhi,absflag);
% cf: filter center frequencies on ERB-rate scale 
% expat: excitation pattern 
% y: waveform vector
% rate: sample rate in Hz (default 10000 Hz)
% nfilt: number of filter channels (default 128)
% flo, fhi: CF limits (default 80 and rate/2=5000 Hz)
% absflag: set to 0 if absolute threshold 
% compensation is NOT desired. (default 1)
% See also: MakeERBFilters, ERBFilterBank

% Peter Assmann - 16-Oct-96
function [cf,expat]=gtep(y,rate,nfilt,flo,fhi,absflag);

npts=max(size(y));
if ~exist('nfilt') nfilt=128; end;
if ~exist('rate') rate=10000; end;
if ~exist('flo') flo=80; end;
if ~exist('fhi') fhi=rate/2; end;
if ~exist('absflag') absflag=1; end;

EarQ = 9.26449;               %  Glasberg and Moore Parameters
minBW = 24.7;
order = 1;
cf = -(EarQ*minBW) + exp((1:nfilt)'*(-log(fhi + EarQ*minBW) + ...
                            log(flo + EarQ*minBW))/nfilt) ...
                          *(fhi + EarQ*minBW);

fcoefs=MakeERBFilters(rate,nfilt,flo);
cf=flipud(cf); % convert CF to ascending order of frequency

% absolute threshold function assumed to adjust the level 
% at the INPUT to the auditory filters 
% (see Moore & Glasberg 1987, p. 218)
% Design a minimum phase filter to apply absthr correction:
if absflag,
 nbin=1025; 
 npts=2*(nbin-1); 
 odd = fix(rem(npts,2)); 
 xx=log(max(sqrt(absthr2(cf)),1e-20));
 xhat = 2*real(ifft(xx,npts));
 wn=[1; 2*ones((npts+odd)/2-1,1); 1; zeros((npts-odd)/2-1,1)];
 absfilt = real(ifft(exp(fft(wn.*xhat(:),npts)),npts));
 y=conv(y,absfilt);
end;

f=flipud(ERBFilterBank(y,fcoefs))';
[npts,nfilt]=size(f);
logarg=sum(f.^2)/npts;
expat=10*log10(max(logarg,eps));
expat=expat+3; % Ad hoc correction factor of 3 dB