% Input is in the frequency domain.  Provide the frequency and amplitude
% of each component in the stimulus.  Save in an m-file as a 2 x N matrix where
% row 1 = the frequencies of the components in Hz, and row 2 = the amplitudes in dB
component = load('MG.m');
component_freq = component(1,:);
component_dB = component(2,:);

% Provide corrections for transduction through outer and middle ear
% and for transducer type (e.g. TDH 39, EAR insert, free field)
% Save corrections in an m-file as a 2 X N matrix where
% row 1 = values in Hz, and row 2 = the corrections in dB SPL
correction = load('MAPC.m');
correction_freq = correction(1,:);
correction_value = correction(2,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PARAMETERS
% c1,c2,c3 used in transforming erb & frequency scales
c1 = 24.673;  c2 = 4.368;  c3 = 2302.6/(c1 * c2);
% range (in erb-rate) over which excitation calculated
eend = 36;  estart = 3;
estep = 0.1; %erb spacing (i.e. separation of adjacent filters)
p51_1k = 4 * 1000/(c1*(c2 + 1)); %for symmetrical filters at 1 kHz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Corrects the level of each component for transducer/transduction characteristics
component_dB = component_dB - interp1(correction_freq, correction_value,...
   component_freq,'cubic');

% Creates an N x M matrix where N = number of component frequencies and
% M =(eend-estart) * estep
erb_steps = ((estart:estep:eend)' * ones(1,size(component_dB,2)))';

% Converts to dB to linear intensity
component_intensity = (10.^(component_dB / 10))' * ones(1,size(erb_steps,2));
clear component_dB

% Converts frequency in Hz to corresponding E-value
component_erb = (c3 * log10(c2 .* (component_freq./1000) + 1))'...
   * ones(1,size(erb_steps,2));

% Computes total intensity in each ERB & converts intensity sum to dB
erb_dB = sum(component_intensity .* (component_erb > erb_steps - (.5 * erb_steps)...
   & component_erb < erb_steps + (.5 * erb_steps)));
erb_dB = ((10 * log10(erb_dB + (erb_dB < exp(-10))* exp(-10)))'...
   * ones(1,size(component_erb,1)))';
clear component_dB

% Computes matrices for the following:
% G (normalized deviation of filter center frequency from each component)
% P51 (used when filter center frequency is below the component--positive G)
% P (used when filter center frequency is above the component--negative G)
freq_Hz = 1000 * (10.^(erb_steps/c3) - 1)/ c2;
component_freq = component_freq' * ones(1,length(erb_steps));
erb = c1 * (c2 * (freq_Hz / 1000) + 1);
P51 = (freq_Hz ./ erb) * 4;
clear erb
P = P51 - .35 .* (P51 / p51_1k) .* (erb_dB - 51);
clear erb_dB
G = (component_freq - freq_Hz) ./ freq_Hz;
clear component_freq

% Calculates the excitation pattern as the total output from
% each auditory filter as a function of filter center frequency (G)
Excitation_Level = 10 * log10(sum(((G > 2) .* exp(-10)) + ((G >= 0 & G <= 2)...
   .* ((1 + P51 .* G) .* exp(-P51 .* G).* component_intensity)) + ((G < 0)...
   .* ((1 + P .* abs(G)) .* exp(-P .* abs(G)).* component_intensity))));
clear component_intensity

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Plot the excitation pattern on the erb-rate scale
figure
plot(erb_steps, Excitation_Level)
AXIS([estart eend 0 (max(Excitation_Level) + 5)])
TITLE('Excitation Pattern')
XLABEL('ERB-rate')
YLABEL('dB Excitation')

% Alternate Plot:
% excitation pattern on a linear frequency scale

%figure
%plot(freq_Hz, Excitation_Level)
%AXIS([freq_Hz(1,1) freq_Hz(1,length(freq_Hz)) 0 (max(Excitation_Level) + 5)])
%TITLE('Excitation Pattern')
%XLABEL('Frequency (Hz)')
%YLABEL('dB Excitation')