function f=gammatones(x,f)
%GAMMATONES Performs gammatone filtering on input signal.
%
%    filterspec=gammatones(x,filterspec)
%
% where filterspec specifies a filterbank (or just an individual
% filter), usually created with function initgammatones.
%
% If x is a single sample, this function performs one sample of
% filtering and returns the result in structure filterspec. If x contains
% a vector, it is all filtered. However, in the latter case, only the 
% most recent values of filtering parameters are maintained. Thus, if you
% need access to these values, either call this function one sample at a time,
% or write a new function.

for t=1:length(x)
  f.tim = f.tim + 1;
  f.p = [ x(t)*cos(f.tim.*f.ec) + sum(f.a.*f.p); f.p(1:3,:)];
  f.q = [-x(t)*sin(f.tim.*f.ec) + sum(f.a.*f.q); f.q(1:3,:)];
end

% now calculate parameters

tmpa = [ones(1,length(f.cf)); f.a(1,:); f.a5];
u0 = sum(tmpa.*f.p(2:4,:));  
v0 = sum(tmpa.*f.q(2:4,:));
f.bm=f.gain.*(u0.*cos(f.tim.*f.ec)-v0.*sin(f.tim.*f.ec));