%%
cd '/Volumes/VOLUME2_MAC_INT_DH/D_DRIVE/CURRENT_WORK/PROMET_2009_MATLAB'
	
% #K
% #G
% DIMENSION_FREQMAP: 30.0;
% DIMENSION_AUDMAP: 4.0;
% DIMENSION_EYEPOS: 10.0;
% DIMENSION_HIDDEN: 10.0;
% DIMENSION_OUTPUT: 20.0;
% 
% SIGMA_FREQMAP: 0.5;
% SIGMA_MOTMAP: 10.0;
% 
% GAIN_FREQMAP: 1;
% EXP_FREQMAP: 40;
% 
% GAIN_EYEPOS: 1;
% GAIN_MOTOR: .7;
% 
% BIAS_HID: 5;
% 
% LEARNING_RATE: 0.5;
% 
% AZIMUTH: 0.0;
% EYEPOSITION: 0.0;
% FREQUENCY: 4950.0;
% INTENSITY: 0.0;
% 
% RANDOM_FLAG: 1;
% CONTINUE_FLAG: 0;
% 
% MOTORERROR_FLAG: 0;
% AZIMUTH_FLAG: 0;
% EYE_FLAG: 0;
% FREQUENCY_FLAG: 0;
% INTENSITY_FLAG: 0;
% 
% DEMO_FLAG: 1000;
% DEMO_WAIT: 100;
% STORE_FLAG: 1;
% STORE: 0;
% 
% WEIGHTS_RANGE: .5;
% 
% AZIMUTH_RANGE: 60.0;
% EYEPOS_RANGE: 40.0;
% MOTORERROR_RANGE: 100.0:
% LOW_FREQUENCY_RANGE: 100.0;
% HIGH_FREQUENCY_RANGE: 10000.0;
% INTENSITY_RANGE: 0.40;
% 
% #V
% 1/1.wht;
% 2/2.wht;
% 3/3.wht;
% 4/4.wht;
% 5/4.wht;
% #E
% 
% 
% 



% 	/* *************** SETTING DIMENSION VALUES ************** */
% 
% 	 Nfreq	 = (unsigned short)get_bat(bt, 0, "DIMENSION_FREQMAP", &Found);
% 	 Naud	 = (unsigned short)get_bat(bt, 0, "DIMENSION_AUDMAP", &Found);
% 	 Nmsc	 = (unsigned short)get_bat(bt, 0, "DIMENSION_EYEPOS", &Found);
% 	 Nmot	 = (unsigned short)get_bat(bt, 0, "DIMENSION_OUTPUT", &Found);
% 	 Nhidden = (unsigned short)get_bat(bt, 0, "DIMENSION_HIDDEN", &Found);
% 
% 
% 	 Sigma	 = get_bat(bt, 0, "SIGMA_MOTMAP", &Found);
% 	 Sig_freq= get_bat(bt, 0, "SIGMA_FREQMAP", &Found);
% 
% 	 Gain_freq= get_bat(bt, 0, "GAIN_FREQMAP", &Found);
% 	 Gain_pos = get_bat(bt, 0, "GAIN_EYEPOS", &Found);
% 	 Gain_mot = get_bat(bt, 0, "GAIN_MOTOR", &Found);
% 	 Exp_freq = get_bat(bt, 0, "EXP_FREQMAP", &Found);
% 
% 	 Eye_range = get_bat(bt, 0, "EYEPOS_RANGE", &Found);
% 	 Mot_range = get_bat(bt, 0, "MOTORERROR_RANGE", &Found);
% 	 Azi_range = get_bat(bt, 0, "AZIMUTH_RANGE", &Found);
% 	 Freq_min  = get_bat(bt, 0, "LOW_FREQUENCY_RANGE", &Found);
% 	 Freq_max  = get_bat(bt, 0, "HIGH_FREQUENCY_RANGE", &Found);
% 	 Int_range = get_bat(bt, 0, "INTENSITY_RANGE", &Found);
% 
% 
% 
% 	 LEARNING_RATE = get_bat(bt, 0, "LEARNING_RATE", &Found);
% 	 BIAS_HID = get_bat(bt, 0, "BIAS_HID", &Found);
% 	 Weights_range = get_bat(bt, 0, "WEIGHTS_RANGE", &Found);
% 
% 	 MOTORERROR_FLAG = (short)get_bat(bt, 0, "MOTORERROR_FLAG", &Found);
% 	 AZIMUTH_FLAG = (short)get_bat(bt, 0, "AZIMUTH_FLAG", &Found);
% 	 EYE_FLAG = (short)get_bat(bt, 0, "EYE_FLAG", &Found);
%  	 FREQUENCY_FLAG = (short)get_bat(bt, 0, "FREQUENCY_FLAG", &Found);
% 	 INTENSITY_FLAG = (short)get_bat(bt, 0, "INTENSITY_FLAG", &Found);
% 	 RAND_FLAG = (short)get_bat(bt, 0, "RANDOM_FLAG", &Found);
% 
% 	 CONTINUE_FLAG = (short)get_bat(bt, 0, "CONTINUE_FLAG", &Found);
% 	 STORE_FLAG = (short)get_bat(bt, 0, "STORE_FLAG", &Found);
% 	 DEMO_FLAG = (short)get_bat(bt, 0, "DEMO_FLAG", &Found);
% 	 STORE = (unsigned long)get_bat(bt, 0, "STORE", &Found);
% 	 WAIT = (unsigned long)get_bat(bt, 0, "DEMO_WAIT", &Found);
% 
% 
% 	/* ************* CREATING VECTORS AND MATRICES *********** */
% 
% 	 Audmap        = vector(1, Naud);
% 	 Eyepos        = matrix(1, 2, 1, Nmsc);
% 	 Cochlea       = matrix(1, 2, 1, Nfreq);
% 
% 	 Hidden        = vector(1, Nhidden);
% 
% 	 MotorError    = vector(1, Nmot);
% 	 Teacher       = vector(1, Nmot);
% 	 Error_hm      = vector(1, Nmot);
% 	 Error_ah      = vector(1, Nhidden);
% 
% 	 Bias_HID      = vector(1, Nhidden);
% 	 Bias_MOT      = vector(1, Nmot);
% 	 Bias_AUD      = vector(1, Naud);
% 
% 	 Weights_AH    = matrix(1, Nhidden, 1, Naud);
% 	 Weights_CA    = matrix(1, Naud, 1, 2 * Nfreq);
% 	 Weights_EH    = matrix(1, Nhidden, 1, 2 * Nmsc);
% 	 Weights_HM    = matrix(1, Nmot, 1, Nhidden);
% 
% /* ************************************************************************************************* */
% 

%%
% DIMENSION_FREQMAP: 30.0;
% DIMENSION_AUDMAP: 4.0;
% DIMENSION_EYEPOS: 10.0;
% DIMENSION_HIDDEN: 10.0;
% DIMENSION_OUTPUT: 20.0;
% 
% SIGMA_FREQMAP: 0.5;
% SIGMA_MOTMAP: 10.0;
% 
% GAIN_FREQMAP: 1;
% EXP_FREQMAP: 40;
% 
% GAIN_EYEPOS: 1;
% GAIN_MOTOR: .7;
% 
% BIAS_HID: 5;
% 
% LEARNING_RATE: 0.5;

		     

%%
% /* ******************************* CALC COCHLEA-LAYER ACTIVATION VALUES ****************** */
% /*
%    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] )
%        */


%void coch_input(cochlea, azi, nstimuli, Is, freq_flag)
clc
freq_flag=0;
Nfreq=30;
Freq_min=100;
Freq_max=10000;

Gain_freq = 1.0;
Exp_freq  = 40.0;
Sig_freq  = 0.5;


azi=30;
nstimuli=4;
Is=1;

    
%%

azi=-90:1:90;
IL=Is-exp( +(1/Exp_freq)*(azi - 90))
IR=Is-exp( -(1/Exp_freq)*(azi + 90))
figure(1)
subplot(211)
plot(azi,IL)
subplot(212)
plot(azi,IR)

%%
clc

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
	 
	 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');
	 pause(0.001)

	
%%
%%%%%% AS FUNCTION OF CF
clc

cf=1000
azi=-60;
Is=1.0;
Gain_freq = 1.0;
Exp_freq  = 40.0;
Sig_freq  = 5;
Nfreq=100;

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
 
 for cf=Freq_min:20:Freq_max
	 io = (log(cf)-log(Freq_min)) / (log(Freq_max)-log(Freq_min)) * (Nfreq -1) +1;
	 cochlea=zeros(2,Nfreq);
	 for i=1:1:Nfreq
		 
		 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
	 
	 plot((cochlea(1,:)),'k');
	 hold 
	 plot((cochlea(2,:)),'r');
	 pause(0.001)
 end
%%


 cochlea=zeros(2,Nfreq);

	 %/* determine cochlear location in frequency map */

	 m = (Nfreq - 1) / (log(Freq_max) - log(Freq_min));

		for k=1:1:nstimuli
				
				if (not(freq_flag))
		 			i0 = randi(Nfreq,1,1);
					m = (Nfreq - 1) / (log(Freq_max) - log(Freq_min));
					logf =  ((i0 - 1.0)/m) + log(Freq_min);
					c_fr =  exp(logf);
				end
					
				if (freq_flag)
					i0=15;
					c_fr = (Freq_max - Freq_min) / 2.0;
				end
					
				
				for i=1:1:Nfreq	
				% 	 /* nonlinear, saturating, azimuth dependent activity for both cochleae:
                % 	    the threshold is put at 90 deg left and right */
                if (azi < -90.0) rate_r = 0.0;
                    else rate_r = Is + (Gain_freq * (1 - exp( -(azi + 90.0) / Exp_freq) ))
				end
				
	            if (azi > 90.0)  rate_l = 0.0;
	                else rate_l = Is + (Gain_freq * (1 - exp( +(azi - 90.0) / Exp_freq) ));
				end	
				
				if(rate_r<0) rate_r=0;end
				if(rate_l<0) rate_l=0;end
					
				arg = ( sqrt(i-i0) / (2 * sqrt(Sig_freq)) );
				
				cochlea(1,i) = cochlea(1,i) + rate_r * exp(-arg)
				if (cochlea(1,i) > rate_r)  cochlea(1,i) = rate_r; end
				
				cochlea(2,i) = cochlea(2,i) + rate_l * exp(-arg)
				if (cochlea(2,i) > rate_l)  cochlea(2,i) = rate_l; end
				end
				
				
		end % /* end k loop */
	%/* *************** END COCHLEA OUTPUT *************** */

	
	
	
	figure(1)
	plot((cochlea(1,:)))
	
	
%%
	
	
	
	
	
	
	
	
	
