% enter with xin = filename
% file is assumed to be a binaural impulse response with the source on the left side
[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
box_length = round(100*sr/1000);  % 100ms window
early_time = round(7*sr/1000); %7ms direct sound
D = box_length;       %the window width
wb = [2*800/sr 2*4000/sr]; % filter the input
[b a] = butter(3,wb);
% wb = [2*1000/sr 2*4000/sr]; % másik szűrő (régi)
% [b a] = ellip(3,2,30,wb); % másik szűrő (régi)
ir_left = filter(b,a,ir_left); ir_right = filter(b,a,ir_right);
datamax = max(ir_left);
for index1 = 1:0.1*sr
   if abs(ir_left(index1))+abs(ir_right(index1)) > datamax/100
       break
   end
end

% for il = 1:0.1*sr % másik (régi) loop
% 	if abs(ir_left(il)) > 500
% 		break
% 	end
% 	if abs(ir_right(il)) > 500
% 		break
% 	end
% end

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

S = upper_scale-10*log10(sum(ir_left.^2));
early = 10*log10(sum(ir_left(1:early_time).^2));
ln = length(ir_left);
log_rvb = 10*log10(cumsum(ir_left(early_time:ln).^2));
for ix = 1:ln-early_time % look at positive values of S+log_rvb
	if S+log_rvb(ix) < 0
		log_rvb(ix) = -S;
	end
end
LOC = -1.5 +(early+S) - (1/D)*sum(S+log_rvb(1:D-early_time));

% the limits look different from the equation – but this is actually OK.
%************************* 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; 
left_power = sqrt(sum(ir_left_rvb.^2));
right_power = sqrt(sum(ir_right_rvb.^2));
n = round(sr*122/1000);
n2 = round(sr*200/1000);
n3 = box_length;
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_left = 20*log10(build_up_left(ln));
direct_reverb_left = 20*log10(direct_level_left)- assymptote_left
log_build_up_left = 20*log10(build_up_left)- assymptote_left;
direct_level_left_plot(1:n) = direct_reverb_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)), axis([0 120 -22 5])
hold on
plot(1000*(1:n2)/sr,log_build_up_left(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('left_channel time in ms')
ylabel('rate of nerve firings - sound energy in dB')
