clear 
fs=48000; % frequency of sampling
%sig=(rand(1,fs/2)-0.5)*10; % 500 ms of white noise
sig=sin([1:fs/2]*1000);

%%
winlen=round(fs/40); % 25 ms time window
[blp,alp]=butter(2,(500 / (fs/2)), 'high'); % high-pass filter
zE=[1:41];                    % utilized ERB channel numbers
fE=228.7*(10.^(zE/ 21.3)-1);  % corresponding frequencies
% gain to implement hump around 4 kHz (ERB 25 +- 7)
hump_coeffs=1+(max(0,7-abs(25-zE))/7)*6;
% approximated spreading function of excitation in ERB scale
spreadfunct=10.^([-80 -60 -40 -20 0 -8 -16 -24 -32 -40 -48 -56  -64]/10);

%%

sig=filter(blp,alp,sig);%high-pass to simulate LF sensitivity loss
a=1; % time position counter

%%
for i=1:winlen/2:(length(sig)-winlen) % loop through the signal
                                                  
 % window the signal, and take FFT
    SIG=fft(hamming(winlen) .* sig(i:(i+winlen-1))');
 
    
    POWSPECT=SIG.*conj(SIG);     % compute power spectrum
    


    % scaling linear frequency to ERB
    lowlimit=1; i=1;
    for z=zE(2:end)
        highlimit=round(fE(z)/fs*winlen); % upper FFT-bin for ERB channel
        % sum the power inside ERB channel
        excitation(i)=sum(POWSPECT(lowlimit:highlimit))*hump_coeffs(i);
        lowlimit=highlimit+1; i=i+1; % update counters and lowlimit
    end


    % implement excitation spreading by convolution with spreading
    % function and store the excitation patter
    excitpattern(a,:)=conv(excitation,spreadfunct);
    

        
    a=a+1; %counter update
end
% computation of auditory spectrum
audspec=10*log10(mean(excitpattern,1)); % avg over time


hearthr=(zE-24).^2/15;  % crude approx. hearing threshold
figure(1); clf;  axes('Position',[0.1 0.1 0.5 0.3])
plot(zE(2:end)-0.5, max(hearthr(2:end), audspec(5:(end-8))),'-')
hold on;  plot(zE(1:end), hearthr,'--'); set(gca,'XTick',[2:4:40]);
xlabel('ERB scale'); ylabel('Auditory spectrum [dB]')