%% Create Gammatone Spectrogram
% Read in an audio signal and convert it to mono for easy visualization.

[audio,fs] = audioread('WaveGuideLoopOne-24-96-stereo-10secs.aif');
audio = mean(audio,2);

%
% Create a gammatoneFilterBank with 64 filters that span the range 62.5 to 20,000 Hz. Pass the audio signal through the filter bank.

gammaFiltBank = gammatoneFilterBank('SampleRate',fs, ...
                                    'NumFilters',64, ...
                                    'FrequencyRange',[62.5,20e3]);

audioOut = gammaFiltBank(audio);

%
% Calculate the energy-per-band using 50 ms windows with 25 ms overlap. Use dsp.AsyncBuffer to divide the signals into overlapped windows and then to log the RMS value of each window for each channel.

samplesPerFrame = round(0.05*fs);
samplesOverlap = round(0.025*fs);

buff = dsp.AsyncBuffer(numel(audio));
write(buff,audioOut.^2);

sink = dsp.AsyncBuffer(numel(audio));

while buff.NumUnreadSamples > 0
    currentFrame = read(buff,samplesPerFrame,samplesOverlap);
    write(sink,mean(currentFrame,1));
end

%
% Convert the energy values to dB. Plot the energy-per-band over time.

gammatoneSpec = read(sink);
D = 20*log10(gammatoneSpec');

timeVector = ((samplesPerFrame-samplesOverlap)/fs)*(0:size(D,2)-1);
cf = getCenterFrequencies(gammaFiltBank)./1e3;

surf(timeVector,cf,D,'EdgeColor','none')
axis([timeVector(1) timeVector(end) cf(1) cf(end)])
view([0 90])
caxis([-150,-60])
colorbar
xlabel('Time (s)')
ylabel('Frequency (kHz)')

%
% Copyright 2018 The MathWorks, Inc.