function [e,time]=auditoryfbank(x,fs,nframeshz,cuts)
% auditoryfbank performs an auditory filterbank analysis on a signal
%   [e,time]=auditoryfbank(x,fs,nframeshz,cuts) 
%   Inputs:   x: signal (trace) 
%             fs: sampling frequency  Hz samples/sec
%             nframeshz: energies to nframeshz  (frames per second)
%                 energies in each bank is output this many times per
%                 second
%             cuts: filter bank frequency cutoffs
%   through a n-channel auditory filterbank, smoothing and downsampling
%   The size of the filter bank is set by the vector cuts
%   Output is in decibels power in each band +100
%   returns an array of e=db,channel and a time vector for the samples
%   Maximum Frequency in cuts must not be larger than fs/2=Nyquist
%   time seems to be in seconds from beginning of file

% get output energy sample points
x=reshape(x,[length(x) 1]);         % force to be a column
%nframeshz = 1;                       % how many frames a second to sample
eidx=int32(fs/nframeshz):int32(fs/nframeshz):length(x); % sampling indices 
ne=length(eidx);                    % number of output energy frames

% these are the filter cut-offs in Hz
% cuts=[0.25 2 5 10 18 25];  % note that maximum must not be larger than fs/2
nf=length(cuts)-1;                  % number of filters

% initialise the output energy table
e=zeros(ne,nf);
time = 0.0:1.0/nframeshz:(length(x)/fs - 1.0/nframeshz);  % time of samplings

% build a smoothing filter
[sb,sa]=butter(4,nframeshz/fs);   % low-pass at half output frame rate
% 4 is order of butterworth filter
% the second argument is cuttoff for lf filter
% 1 is the nyquist frequency which would be 0.5fs

% for each channel in turn, nf is number of filters
for i=1:nf
    fprintf('Calculating channel %d (%d-%dHz)\n',i,cuts(i),cuts(i+1));
    % build the band-pass filter
    [b,a]=butter(4,[2*cuts(i)/fs 0.999999*2*cuts(i+1)/fs]);
    % filter the signal this filter
    y=filter(b,a,x);
    % rectify and smooth with the smoothing filter
    sy=filter(sb,sa,abs(y));
    % sample and convert to decibels
    sy=max(sy,eps);                   % force to be > 0, eps is smallest no
    e(:,i)=100+20*log10(sy(eidx));    % downsample and convert to dB
end

 