%
%     gammatone
%     gt(t) = t^(n-1)*exp(-2pi*bt)*cos(2pi*ft)
%     parameters: n, b, and f
%     (default: n=4, b=125 Hz and f = 1000 Hz)
%  
function gammatone(n, b, cf)
if nargin<3
    cf=1000;
end
if nargin<2
    % b=125;
    b=24.7*(1+4.37*cf/1000);   % Cambridge critical bandwidth
end
if nargin<1
    n=4;
end

t=0:0.00001:0.035;
gt = t.^(n-1).*exp(-2*pi*b*t).*cos(2*pi*cf*t);

figure(1)
clf
plot(t,gt,'k-','linewidth',1.5);
hold on
set(gca, 'fontsize',12,'fontweight','bold');
xlabel('time (s)', 'fontsize',13, 'fontweight','bold');
ylabel('amplitude (a.u.)', 'fontsize',13, 'fontweight','bold');
mx=max(gt);
s=sprintf('n=%1d  b=%3.0f Hz   CF = %4d Hz',n,b,cf);
text(0.0075, 0.9*mx, s, 'fontsize',12,'fontweight','bold');
axis([0 0.015 -1.05*mx 1.05*mx]);

print -djpeg gammatone

%
%   the amplitude spectrum of the filter
%
fg = fft(gt, 4096);
fga=abs(fg(1:2048));
mxf = max(fga);
fga=fga/mxf;
FMAX = 0.5/0.00001;
freqs = 0:FMAX/2048:FMAX-FMAX/2048;
figure(2)
clf
loglog(freqs, fga, 'k-','linewidth',2);
hold on
set(gca, 'YDIR','reverse');
set(gca, 'fontsize',12,'fontweight','bold');
xlabel('Frequency (Hz)', 'fontsize',13,'fontweight','bold');
ylabel('Amplitude re. maximum', 'fontsize',13,'fontweight','bold');
axis([0 20000 10^(-6) 1.2])
%s=sprintf('Gamma filter: n=%1d  b=%3.0f Hz   CF = %4d Hz',n,b,cf);
%text(30,10^(-6),s,'fontsize',12,'fontweight','bold');
%plot([500 2000], [10^(-1) 10^(-1)], 'k--','linewidth',1.5);

%
%   the phase spectrum of the filter
%
fgf1=angle(fg(1:2048));
fgf = unwrap(fgf1);
figure(3); clf;
semilogx(freqs, fgf/pi+1,'k-','linewidth',2);
hold on
set(gca, 'fontsize',12,'fontweight','bold');
xlabel('Frequency (Hz)', 'fontsize',13,'fontweight','bold');
ylabel('Phase (in pi rad)', 'fontsize',13,'fontweight','bold');
axis([0 20000  -3 1])

for f=[125, 250, 500, 2000, 4000]
    b=24.7*(1+4.37*f/1000);
    gt = t.^(n-1).*exp(-2*pi*b*t).*cos(2*pi*f*t);
    fg = fft(gt, 4096);
    fga=abs(fg(1:2048));
    mxf = max(fga);
    fga=fga/mxf;
    figure(2);
    loglog(freqs, fga, 'k--','linewidth',1.0);
end

print -djpeg gammatonespectrum

for f=[125, 250, 500, 2000, 4000]
    b=24.7*(1+4.37*f/1000);
    gt = t.^(n-1).*exp(-2*pi*b*t).*cos(2*pi*f*t);
    fg = fft(gt, 4096);
    fgf1=angle(fg(1:2048));
    fgf=unwrap(fgf1);
    figure(3);
    semilogx(freqs, fgf/pi+1, 'k--','linewidth',1.0);
end

print -djpeg gammatonephasespectrum


%
%    fc = 500 Hz
%
cf = 500;
b=24.7*(1+4.37*cf/1000); 

t=0:0.00001:0.035;
gt = t.^(n-1).*exp(-2*pi*b*t).*cos(2*pi*cf*t);

figure(4)
clf
plot(t,gt,'k-','linewidth',1.5);
hold on
set(gca, 'fontsize',12,'fontweight','bold');
xlabel('time (s)', 'fontsize',13, 'fontweight','bold');
ylabel('amplitude (a.u.)', 'fontsize',13, 'fontweight','bold');
mx=max(gt);
s=sprintf('n=%1d  b=%3.0f Hz   CF = %4d Hz',n,b,cf);
text(0.0075, 0.9*mx, s, 'fontsize',12,'fontweight','bold');
axis([0 0.025 -1.05*mx 1.05*mx]);

print -djpeg gammatone500

%
%   the amplitude spectrum of the filter
%
fg = fft(gt, 4096);
fga=abs(fg(1:2048));
mxf = max(fga);
fga=fga/mxf;
figure(5)
clf
subplot(1,2,1);
loglog(freqs(1:60), fga(1:60), 'k-','linewidth',2);
set(gca, 'YDIR','reverse');
hold on
set(gca, 'fontsize',12,'fontweight','bold');
xlabel('Frequency (Hz)', 'fontsize',13,'fontweight','bold');
ylabel('Amplitude re. maximum', 'fontsize',13,'fontweight','bold');
axis([0 1100 10^(-4.5) 1.2])
set(gca,'xtick',[0 100 200 300 400 500 600 700 800 900 1000 1100]);
set(gca, 'xticklabel', ['    ';' 100';'    ';'    ';'    ';' 500';'    ';'    ';'    ';'    ';'1000']);

fgf1=angle(fg);
fgf = unwrap(fgf1);
subplot(1,2,2)
semilogx(freqs(1:60), fgf(1:60)/pi+1,'k-','linewidth',2);
hold on
set(gca, 'fontsize',12,'fontweight','bold');
xlabel('Frequency (Hz)', 'fontsize',13,'fontweight','bold');
ylabel('Phase (in pi rad)', 'fontsize',13,'fontweight','bold');
axis([0 1100  -3 1])
set(gca,'xtick',[0 100 200 300 400 500 600 700 800 900 1000 1100]);
set(gca, 'xticklabel', ['    ';' 100';'    ';'    ';'    ';' 500';'    ';'    ';'    ';'    ';'1000']);

print -djpeg gammatonespectrum500

%
%   the phase spectrum of the filter
%
fgf1=angle(fg);
fgf = unwrap(fgf1);
figure(6); clf;
semilogx(freqs(1:60), fgf(1:60)/pi+1,'k-','linewidth',2);
hold on
set(gca, 'fontsize',12,'fontweight','bold');
xlabel('Frequency (Hz)', 'fontsize',13,'fontweight','bold');
ylabel('Phase (in pi rad)', 'fontsize',13,'fontweight','bold');
axis([0 10000  -3 1])

print -djpeg gammatonephasespectrum500

return
