function bm=gammatone(x,fshz,cfhz)

%GAMMATONE gammatone(x,fshz,cfhz,pc) implements the 4th 
% order gammatone filter with centre frequency cfhz (in Hz)
% and sampling frequency fshz (in Hz).
% If phase-correction is required, supply anything as the
% 4th argument.
% If called with zero or one output arguments, it returns the 
% filter output i.e. simulated basilar membrane displacement. 
%
% Called as [bm,env] = gammatone(x,fshz,cfhz) it returns the 
% filter output and instantaneous envelope. 
%
% Called as [bm,env,instf] = gammatone(x,fshz,cfhz) it returns 
% the filter output, instantaneous envelope and frequency.
%

% Implementation by impulse invariant transform as described in
% Cooke's thesis.
% Extended by Martin from Guy's version - 20 Feb 97
% zero-pad if phase-correction is required
% cut down from Martins version by guy!

wcf=2*pi*cfhz;									% radian frequency
tpt=(2*pi)/fshz;
bw=erb(cfhz)*bwcorrection;
a=exp(-bw*tpt);
kT=[0:length(x)-1]/fshz;
gain=((bw*tpt)^4)/6;							% based on integral of impulse response
q=exp(-j*wcf*kT).*x;							% shift down to d.c.
p=filter([1 0],[1 -4*a 6*a^2 -4*a^3 a^4],q);	% filter: part 1
u=filter([1 4*a 4*a^2 0],[1 0],p); 				% filter: part 2
bm=gain*real(exp(j*wcf*kT).*u);					% shift up in frequency