function [Y,F] = gammatonegram(X,SR,TWIN,THOP,N,FMIN,FMAX,USEFFT,WIDTH)
% [Y,F] = gammatonegram(X,SR,N,TWIN,THOP,FMIN,FMAX,USEFFT,WIDTH)
% Calculate a spectrogram-like time frequency magnitude array
% based on Gammatone subband filters. Waveform X (at sample
% rate SR) is passed through an N (default 64) channel gammatone
% auditory model filterbank, with lowest frequency FMIN (50)
% and highest frequency FMAX (SR/2). The outputs of each band
% then have their energy integrated over windows of TWIN secs
% (0.025), advancing by THOP secs (0.010) for successive
% columns. These magnitudes are returned as an N-row
% nonnegative real matrix, Y.
% If USEFFT is present and zero, revert to actual filtering and
% summing energy within windows.
% WIDTH (default 1.0) is how to scale bandwidth of filters
% relative to ERB default (for fast method only).
% F returns the center frequencies in Hz of each row of Y
% (uniformly spaced on a Bark scale).
%
% 2009-02-18 DAn Ellis dpwe@ee.columbia.edu
% Last updated: $Date: 2009/02/23 21:07:09 $
if nargin < 2; SR = 16000; end
if nargin < 3; TWIN = 0.025; end
if nargin < 4; THOP = 0.010; end
if nargin < 5; N = 64; end
if nargin < 6; FMIN = 50; end
if nargin < 7; FMAX = SR/2; end
if nargin < 8; USEFFT = 1; end
if nargin < 9; WIDTH = 1.0; end
if USEFFT == 0
% Use malcolm's function to filter into subbands
%%%% IGNORES FMAX! *****
[fcoefs,F] = MakeERBFilters(SR, N, FMIN);
fcoefs = flipud(fcoefs);
XF = ERBFilterBank(X,fcoefs);
nwin = round(TWIN*SR);
% Always use rectangular window for now
% if USEHANN == 1
window = hann(nwin)';
% else
% window = ones(1,nwin);
% end
% window = window/sum(window);
% XE = [zeros(N,round(nwin/2)),XF.^2,zeros(N,round(nwin/2))];
XE = [XF.^2];
hopsamps = round(THOP*SR);
ncols = 1 + floor((size(XE,2)-nwin)/hopsamps);
Y = zeros(N,ncols);
% winmx = repmat(window,N,1);
for i = 1:ncols
% Y(:,i) = sqrt(sum(winmx.*XE(:,(i-1)*hopsamps + [1:nwin]),2));
Y(:,i) = sqrt(mean(XE(:,(i-1)*hopsamps + [1:nwin]),2));
end
else
% USEFFT version
% How long a window to use relative to the integration window requested
winext = 1;
twinmod = winext * TWIN;
% first spectrogram
nfft = 2^(ceil(log(2*twinmod*SR)/log(2)));
nhop = round(THOP*SR);
nwin = round(twinmod*SR);
[gtm,F] = fft2gammatonemx(nfft, SR, N, WIDTH, FMIN, FMAX, nfft/2+1);
% perform FFT and weighting in amplitude domain
Y = 1/nfft*gtm*abs(specgram(X,nfft,SR,nwin,nwin-nhop));
% or the power domain? doesn't match nearly as well
%Y = 1/nfft*sqrt(gtm*abs(specgram(X,nfft,SR,nwin,nwin-nhop).^2));
end