%%
clc



% /* nonlinear, saturating, azimuth dependent activity for both cochleae:
% 	    the threshold is put at 90 deg left and right */
% /*
% This function calculates the activity values of the cochlea layers (left and right) which are presented in a forward
% direction to the audmap layer
% */
% 
% /*
% function that generates Gaussian receptive field profile:
% F = exp( -[(x-x0)^2]/[2sigma^2] )
% */

cf=1000
azi=-60;
Is=1.0;
Gain_freq = 1.0;
Exp_freq  = 40.0;
Sig_freq  = 4;
Nfreq=100;
nstimuli=10;

freq_flag=0;
Nfreq=100;
Freq_min=100;
Freq_max=10000;

% IL=Is-exp( +(1/Exp_freq)*(azi - 90));
% if(IL<0) IL=0.0; end;
% IR=Is-exp( -(1/Exp_freq)*(azi + 90));
% if(IR<0) IR=0.0; end;

 figure(1)
 hold off
 
 cochlea=zeros(2,Nfreq);
 for k=1:1:nstimuli
	 
	 %cf = randi(10000,1,1)+99;
	 %io = (log(cf)-log(Freq_min)) / (log(Freq_max)-log(Freq_min)) * (Nfreq -1) +1
	 
	 %/* determine cochlear location in frequency map */
	 io = randi(Nfreq,1,1)
	 m = (Nfreq - 1) / (log(Freq_max) - log(Freq_min));
	 logf =  ((io - 1.0)/m) + log(Freq_min);
	 cf =  exp(logf);
	 
	 
	 for i=1:1:Nfreq
		 % changed dependency blocking is freqency dependent        %
		 % High freq much diminution                                %
		 % Low frequency less diminution    ((float)Exp_freq+i*20)  %
		 IL=Is-exp( +(1/(Exp_freq + i*20))*(azi - 90));
		 if(IL<0) IL=0.0; end;
		 IR=Is-exp( -(1/(Exp_freq + 1*20))*(azi + 90));
		 if(IR<0) IR=0.0; end;
		 
		 rate_r = Gain_freq  * IR * exp( -((i-io)*(i-io)) / (2*(Sig_freq*Sig_freq)) ); 
		 if (rate_r<IR) rate_L=IR; end;
		 cochlea(1,i) =  cochlea(1,i)+ rate_r;
		 
		 rate_l = Gain_freq  * IL * exp( -((i-io)*(i-io)) / (2*(Sig_freq*Sig_freq)) );
		 if (rate_l<IL) rate_L=IL; end;
		 cochlea(2,i) = cochlea(2,i) + rate_l;
		 	 
	 end
	 
end 
	 plot((cochlea(1,:)),'k');
	 hold on
	 plot((cochlea(2,:)),'r');