% enter with xin = filename
% file is assumed to be a binaural impulse response with the source on the left side
% fudge factors removed from this version. I saw no good reason to keep them.

[Y,sr,bits] = wavread(xin);

ln = length(Y);

ir_left = Y(1:ln,1);
ir_right = Y(1:ln,2);

upper_scale = 20;  % 20dB range for firings

% proposed box length
box_length = round(100*sr/1000);  % try 100ms
early_time = round(5*sr/1000);

D = box_length;       %the window width

wb = [2*800/sr 2*4000/sr];
[b a] = butter(2,wb);

ir_left = filter(b,a,ir_left);
ir_right = filter(b,a,ir_right);

datamax = max(ir_left);

for index1 = 1:0.2*sr
	if abs(ir_left(index1))+abs(ir_right(index1)) > datamax/20
		break
	end
end

ir_left(1:index1-20) = [];
ir_right(1:index1-20) = [];

% ir_left is an ipselateral binaural impulse response, %truncated to start at zero and filtered to 1000-4000Hz.
% early_time is 5ms in samples, D is 100ms in samples.

% here starts the equation on the slide:

S = upper_scale-10*log10(sum(ir_left.^2));

early = 10*log10(sum(ir_left(1:early_time).^2));

% first integral is a cumsum representing the build up in energy when the IR is excited by a steady tone:

ln = length(ir_left);

log_rvb = 10*log10(cumsum(ir_left(early_time:ln).^2));

% look at positive values of S+log_rvb  only

for ix = 1:ln-early_time
	if S+log_rvb(ix) < 0
		log_rvb(ix) = -S;
	end
end
% the log_rvb(1:D-early_time) term is correct, although the equation has different
% limits, because in MATLAB log_rvb starts at early time, not zero

LOC = (early+S) - (1/D)*sum(S+log_rvb(1:D-early_time))

%************************* graph box ********************

ir_left_rvb = ir_left;
ir_left_rvb(1:early_time) = datamax/100000; %zeros(size(1:early_time);
ir_right_rvb = ir_right;
ir_right_rvb(1:early_time) = datamax/100000; %zeros(size(1:early_time);

% this now starts at time zero, and is essentially zero for early_time
left_power = sqrt(sum(ir_left.^2));
right_power = sqrt(sum(ir_right.^2));
total_power = sqrt(sum((ir_left.^2)+sum(ir_right.^2))/2);

%ir_left_rvb = ir_left_rvb/total_power;
%ir_left = ir_left/total_power;
%ir_right_rvb = ir_right_rvb/total_power;
%ir_right = ir_right/total_power;

build_up_total = sqrt(cumsum((ir_left_rvb/total_power).^2+(ir_right_rvb/total_power).^2));

direct_level_total = sqrt(sum((ir_left(1:early_time)/total_power).^2+(ir_right(1:early_time)/total_power).^2));

ln = length(build_up_total);
assymptote_total = 20*log10(build_up_total(ln));

direct_reverb_total = 20*log10(direct_level_total)- assymptote_total
log_build_up_total = 20*log10(build_up_total)- assymptote_total;

n = round(sr*122/1000);
n2 = round(sr*200/1000);
n3 = box_length;

direct_level_total_plot(1:n) = direct_reverb_total;
direct_level_total_plot(n+1:n2) = -80;

window_plot(1:n3) = 0;
window_plot(n3+1:n2) = -80;

zero_line_plot(1:n2) = -upper_scale;

%plot(1000*(1:n2)/sr,direct_level_total_plot(1:n2)), axis([0 120 -22 5])
%hold on
%plot(1000*(1:n2)/sr,log_build_up_total(1:n2),'r')
%plot(1000*(1:n2)/sr,window_plot(1:n2),'k')
%plot(1000*(1:n2)/sr,zero_line_plot(1:n2),'g')
%hold off
%xlabel('both channels time in ms')
%ylabel('rate of nerve firings - sound energy in dB')

%figure

build_up_left = sqrt(cumsum((ir_left_rvb/left_power).^2));
direct_level_left = sqrt(sum((ir_left(1:early_time)/left_power).^2));

ln = length(build_up_left);
assymptote_reverb_left = 20*log10(build_up_left(ln));

assymptote_total_left = 20*log10(left_power);

direct_reverb_left = 20*log10(direct_level_left)- assymptote_reverb_left
log_build_up_left = 20*log10(build_up_left);  %- assymptote_total_left;


direct_level_left_plot(1:n) = 20*log10(direct_level_left);
direct_level_left_plot(n+1:n2) = -80;

window_plot(1:n3) = 0;
window_plot(n3+1:n2) = -80;

zero_line_plot(1:n2) = -upper_scale;

plot(1000*(1:n2)/sr,direct_level_left_plot(1:n2)+20), axis([0 120 0 25])
hold on
plot(1000*(1:n2)/sr,log_build_up_left(1:n2)+20,'r')
plot(1000*(1:n2)/sr,window_plot(1:n2)+20,'k')
%plot(1000*(1:n2)/sr,zero_line_plot(1:n2)+20,'g')
hold off
xlabel('left channel - time in ms')
ylabel('rate of nerve firings - (sound energy in dB)')
title('blue = direct sound    red = build up of reflections')