changeset 34:e097e9044ef6

new compression and ne prob refractory
author Ray Meddis <rmeddis@essex.ac.uk>
date Fri, 19 Aug 2011 16:07:07 +0100
parents 161913b595ae
children 25d53244d5c8
files MAP/MAP1_14.m multithreshold 1.46/MTprofile.m multithreshold 1.46/MTprofile10_0hr19_Aug_2011.m multithreshold 1.46/MTprofile17_14hr18_Aug_2011.m multithreshold 1.46/MTprofile17_32hr18_Aug_2011.m multithreshold 1.46/MTprofile17_33hr18_Aug_2011.m multithreshold 1.46/MTprofile17_44hr18_Aug_2011.m multithreshold 1.46/MTprofile7_44hr19_Aug_2011.m multithreshold 1.46/MTprofile8_19hr19_Aug_2011.m multithreshold 1.46/MTprofile8_24hr19_Aug_2011.m multithreshold 1.46/expGUI_MT.m multithreshold 1.46/paradigms/paradigm_profile.m multithreshold 1.46/plotProfile.m multithreshold 1.46/profile2mFile.m multithreshold 1.46/savedData/11-Aug-2011 10_53_52.mat multithreshold 1.46/savedData/11-Aug-2011 10_55_36.mat multithreshold 1.46/savedData/11-Aug-2011 15_51_34.mat multithreshold 1.46/savedData/11-Aug-2011 15_54_08.mat multithreshold 1.46/savedData/11-Aug-2011 15_57_29.mat multithreshold 1.46/savedData/11-Aug-2011 16_14_22.mat multithreshold 1.46/savedData/11-Aug-2011 16_26_04.mat multithreshold 1.46/savedData/11-Aug-2011 16_47_28.mat multithreshold 1.46/savedData/11-Aug-2011 16_48_53.mat multithreshold 1.46/savedData/11-Aug-2011 16_49_46.mat multithreshold 1.46/savedData/11-Aug-2011 16_50_10.mat multithreshold 1.46/savedData/11-Aug-2011 16_50_37.mat multithreshold 1.46/savedData/12-Aug-2011 05_33_35.mat multithreshold 1.46/savedData/12-Aug-2011 05_41_36.mat multithreshold 1.46/savedData/12-Aug-2011 05_48_34.mat multithreshold 1.46/savedData/12-Aug-2011 05_56_54.mat multithreshold 1.46/savedData/12-Aug-2011 06_40_17.mat multithreshold 1.46/savedData/12-Aug-2011 06_56_59.mat multithreshold 1.46/savedData/12-Aug-2011 07_20_14.mat multithreshold 1.46/savedData/12-Aug-2011 07_29_20.mat multithreshold 1.46/savedData/12-Aug-2011 07_43_17.mat multithreshold 1.46/savedData/12-Aug-2011 07_49_00.mat multithreshold 1.46/savedData/mostRecentResults.mat multithreshold 1.46/temp.m parameterStore/MAPparamsNormal.m profiles/MTprofileMOCis0075.m profiles/MTprofileMOCis01.m testPrograms/testAN.m testPrograms/testANprob.m testPrograms/testBM.m testPrograms/testFM.m testPrograms/testPhysiology.m testPrograms/testPhysiologyProb.m testPrograms/test_MAP1_14.m utilities/Util_timeStamp.m
diffstat 49 files changed, 546 insertions(+), 181 deletions(-) [+]
line wrap: on
line diff
--- a/MAP/MAP1_14.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/MAP/MAP1_14.m	Fri Aug 19 16:07:07 2011 +0100
@@ -204,13 +204,13 @@
 MOCattSegment=zeros(nBFs,reducedSegmentLength);
 MOCattenuation=ones(nBFs,signalLength);
 
-if DRNLParams.a>0
-    DRNLcompressionThreshold=10^((1/(1-DRNLParams.c))* ...
-    log10(DRNLParams.b/DRNLParams.a));
-else
-    DRNLcompressionThreshold=inf;
-end
-
+% if DRNLParams.a>0
+%     DRNLcompressionThreshold=10^((1/(1-DRNLParams.c))* ...
+%     log10(DRNLParams.b/DRNLParams.a));
+% else
+%     DRNLcompressionThreshold=inf;
+% end
+DRNLcompressionThreshold=DRNLParams.cTh;
 DRNLlinearOrder= DRNLParams.linOrder;
 DRNLnonlinearOrder= DRNLParams.nonlinOrder;
 
@@ -380,6 +380,7 @@
 ANprobability=zeros(nANchannels,segmentLength);
 ANprobRateOutput=zeros(nANchannels,signalLength);
 lengthAbsRefractoryP= round(AN_refractory_period/dt);
+cumANnotFireProb=ones(nANchannels,signalLength);
 % special variables for monitoring synaptic cleft (specialists only)
 savePavailableSeg=zeros(nANchannels,segmentLength);
 savePavailable=zeros(nANchannels,signalLength);
@@ -423,6 +424,8 @@
 nCNtauGk=length(CNtauGk);
 
 % the total number of 'channels' is now greater
+% 'channel' is defined as collections of units with the same parameters
+%  i.e. same BF, same ANtau, same CNtauGk
 nCNchannels=nANchannels*nCNtauGk;
 
 nCNneuronsPerChannel=MacGregorMultiParams.nNeuronsPerBF;
@@ -611,9 +614,11 @@
 
         %            *linear* path
         linOutput = stapesDisplacement * linGAIN;  % linear gain
+       
         for order = 1 : GTlinOrder
             [linOutput GTlinBdry{BFno,order}] = ...
-                filter(GTlin_b(BFno,:), GTlin_a(BFno,:), linOutput, GTlinBdry{BFno,order});
+                filter(GTlin_b(BFno,:), GTlin_a(BFno,:), linOutput, ...
+                GTlinBdry{BFno,order});
         end
 
         %           *nonLinear* path
@@ -633,15 +638,48 @@
                 filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ...
                 nonlinOutput, GTnonlinBdry1{BFno,order});
         end
-        %       broken stick instantaneous compression
-        y= nonlinOutput.* DRNLa;  % linear section.
+        
+%         %    original   broken stick instantaneous compression
+%         y= nonlinOutput.* DRNLa;  % linear section.
+%         % compress parts of the signal above the compression threshold
+%         abs_x = abs(nonlinOutput);
+%         idx=find(abs_x>DRNLcompressionThreshold);
+%         if ~isempty(idx)>0
+%             y(idx)=sign(y(idx)).* (DRNLb*abs_x(idx).^DRNLc);
+%         end
+%         nonlinOutput=y;
+
+        
+        %   new broken stick instantaneous compression
+        y= nonlinOutput.* DRNLa;  % linear section attenuation/gain.
         % compress parts of the signal above the compression threshold
-        abs_x = abs(nonlinOutput);
-        idx=find(abs_x>DRNLcompressionThreshold);
+%         holdY=y;
+        abs_y = abs(y);
+        idx=find(abs_y>DRNLcompressionThreshold);
         if ~isempty(idx)>0
-            y(idx)=sign(y(idx)).* (DRNLb*abs_x(idx).^DRNLc);
+%             y(idx)=sign(y(idx)).* (DRNLcompressionThreshold + ...
+%                 (abs_y(idx)-DRNLcompressionThreshold).^DRNLc);
+            y(idx)=sign(y(idx)).* (DRNLcompressionThreshold + ...
+                (abs_y(idx)-DRNLcompressionThreshold)*DRNLc);
         end
         nonlinOutput=y;
+       
+% % Boltzmann compression
+% y=(nonlinOutput * DRNLa*100000);
+% holdY=y;
+% y=abs(y);
+% s=10; u=0.0;
+% x=1./(1+exp(-(y-u)/s))-0.5;
+% nonlinOutput=sign(nonlinOutput).*x/10000;
+
+
+%         if segmentStartPTR==10*segmentLength+1
+%             figure(90)
+%         plot(holdY,'b'), hold on
+%         plot(nonlinOutput, 'r'), hold off
+%         ylim([-1e-5 1e-5])
+%         pause(1)
+%         end
 
         %       second filter removes distortion products
         for order = 1 : GTnonlinOrder
@@ -652,6 +690,7 @@
 
         %  combine the two paths to give the DRNL displacement
         DRNLresponse(BFno,:)=linOutput+nonlinOutput;
+%         disp(num2str(max(linOutput)))
     end % BF
 
     % segment debugging plots
@@ -726,7 +765,7 @@
         CaCurrent=CaCurrent +  ICa(:,idx)*dt - CaCurrent*dt./tauCa;
         synapticCa(:,idx)=CaCurrent;
     end
-    synapticCa=-synapticCa; % treat IHCpreSynapseParams as positive substance
+    synapticCa=-synapticCa; % treat synapticCa as positive substance
 
     % NB vesicleReleaseRate is /s and is independent of dt
     vesicleReleaseRate = synapse_z * synapticCa.^synapse_power; % rate
@@ -759,7 +798,9 @@
                 Preprocess= Preprocess + reuptake - Preprocessed;
                 Pavailable(Pavailable<0)=0;
                 savePavailableSeg(:,t)=Pavailable;    % synapse tracking
+                
             end
+
             % and save it as *rate*
             ANrate=ANprobability/dt;
             ANprobRateOutput(:, segmentStartPTR:...
@@ -767,8 +808,43 @@
             % monitor synapse contents (only sometimes used)
             savePavailable(:, segmentStartPTR:segmentStartPTR+segmentLength-1)=...
                 savePavailableSeg;
+            
+            %% Apply refractory effect
+               % the probability of a spike's occurring in the preceding
+               %  refractory window (t= tnow-refractory period to tnow-)
+               %    pFired= 1 - II(1-p(t)),
+               % we need a running account of cumProb=II(1-p(t))
+               %   cumProb(t)= cumProb(t-1)*(1-p(t))/(1-p(t-refracPeriod))
+               %   cumProb(0)=0
+               %   pFired(t)= 1-cumProb(t)
+               % This gives the fraction of firing events that must be 
+               %  discounted because of a firing event in the refractory
+               %  period
+               %   p(t)= ANprobOutput(t) * pFired(t)
+               % where ANprobOutput is the uncorrected firing probability
+               %  based on vesicle release rate
+               % NB this covers only the absoute refractory period
+               % not the relative refractory period. To approximate this it
+               % is necessary to extend the refractory period by 50%
 
-            % Estimate efferent effects. ARattenuation (0 <> 1)
+
+                for t = segmentStartPTR:segmentEndPTR;
+                    if t>1
+                    ANprobRateOutput(:,t)= ANprobRateOutput(:,t)...
+                        .* cumANnotFireProb(:,t-1);
+                    end
+                    % add recent and remove distant probabilities
+                    refrac=round(lengthAbsRefractoryP * 1.5);
+                    if t>refrac
+                        cumANnotFireProb(:,t)= cumANnotFireProb(:,t-1)...
+                            .*(1-ANprobRateOutput(:,t)*dt)...
+                            ./(1-ANprobRateOutput(:,t-refrac)*dt);
+                    end
+                end
+%                 figure(88), plot(cumANnotFireProb'), title('cumNotFire')
+%                 figure(89), plot(ANprobRateOutput'), title('ANprobRateOutput')
+
+            %% Estimate efferent effects. ARattenuation (0 <> 1)
             %  acoustic reflex
             [r c]=size(ANrate);
             if r>nBFs % Only if LSR fibers are computed
@@ -783,8 +859,7 @@
             end
             %             plot(ARattenuation)
 
-            % MOC attenuation
-            % within-channel HSR response only
+            % MOC attenuation based on within-channel HSR fiber activity
             HSRbegins=nBFs*(nANfiberTypes-1)+1;
             rates=ANrate(HSRbegins:end,:);
             if rateToAttenuationFactorProb<0
@@ -797,9 +872,11 @@
                         filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ...
                         MOCprobBoundary{idx});
                     smoothedRates=smoothedRates-MOCrateThresholdProb;
-                    smoothedRates(smoothedRates<0)=0;
-                    MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= ...
-                        (1- smoothedRates* rateToAttenuationFactorProb);
+                    smoothedRates=max(smoothedRates, 0);
+                    
+                    x=(1- smoothedRates* rateToAttenuationFactorProb);
+                    x=max(x, 10^(-30/20));
+                    MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= x;
                 end
             end
             MOCattenuation(MOCattenuation<0)=0.001;
@@ -1010,11 +1087,12 @@
             %% IC ----------------------------------------------
                 %  MacGregor or some other second order neurons
 
-                % combine CN neurons in same channel, 
-                %  i.e. same BF & same tauCa
+                % combine CN neurons in same channel (BF x AN tau x CNtau) 
+                %  i.e. same BF, same tauCa, same CNtau
                 %  to generate inputs to single IC unit
                 channelNo=0;
-                for idx=1:nCNneuronsPerChannel:nCNneurons-nCNneuronsPerChannel+1;
+                for idx=1:nCNneuronsPerChannel: ...
+                        nCNneurons-nCNneuronsPerChannel+1;
                     channelNo=channelNo+1;
                     CN_PSTH(channelNo,:)=...
                         sum(CN_spikes(idx:idx+nCNneuronsPerChannel-1,:));
@@ -1034,7 +1112,7 @@
                 ICtrailingAlphas=ICcurrentTemp(:, reducedSegmentLength+1:end);
 
                 if IC_c==0
-                    % faster computation when threshold is stable (C==0)
+                    % faster computation when threshold is stable (c==0)
                     for t=1:reducedSegmentLength
                         s=IC_E>IC_Th0;
                         dE = (-IC_E/IC_tauM + inputCurrent(:,t)/IC_cap +...
@@ -1089,8 +1167,13 @@
                 ICoutput(:,reducedSegmentPTR:shorterSegmentEndPTR)=ICspikes;
                 
                 % store membrane output on original dt scale
-                if nICcells==1  % single channel
-                    x= repmat(ICmembranePotential(1,:), ANspeedUpFactor,1);
+                % do this for single channel models only 
+                % and only for the HSR-driven IC cells
+                if round(nICcells/length(ANtauCas))==1  % single channel
+                    % select HSR driven cells
+                    x= ICmembranePotential(length(ANtauCas),:);
+                    % restore original dt
+                    x= repmat(x, ANspeedUpFactor,1);
                     x= reshape(x,1,segmentLength);
                     if nANfiberTypes>1  % save HSR and LSR
                         y=repmat(ICmembranePotential(end,:),...
@@ -1154,6 +1237,18 @@
 
 %% apply refractory correction to spike probabilities
 
+% the probability of a spike's having occurred in the preceding 
+%  refractory window 
+%    pFired= 1 - II(1-p(t)), t= tnow-refractory period: tnow-1
+% we need a running account of cumProb=II(1-p(t))
+%   cumProb(t)= cumProb(t-1)*(1-p(t-1))/(1-p(t-refracPeriod))
+%   cumProb(0)=0
+%   pFired(t)= 1-cumProb(t)
+% whence
+%   p(t)= ANprobOutput(t) * pFired(t)
+% where ANprobOutput is the uncorrected probability
+
+
 % switch AN_spikesOrProbability
 %     case 'probability'
 %         ANprobOutput=ANprobRateOutput*dt;
--- a/multithreshold 1.46/MTprofile.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/multithreshold 1.46/MTprofile.m	Fri Aug 19 16:07:07 2011 +0100
@@ -1,24 +1,31 @@
 function x = MTprofile
+%created: 17_20hr18_Aug_2011
+
 x.BFs = [250   500  1000  2000  4000  8000];
-x.LongTone = [22      17.4      15.7      15.1      20.1        26];
-x.ShortTone = [25.3      22.5      18.4      19.2      25.8      31.8];
+
+x.LongTone = [16.4      13.1      8.77      9.38      15.8      21.4];
+x.ShortTone = [18.8      15.1      11.9      13.3      20.7        25];
+
 x.Gaps = [0.01      0.03      0.05      0.07      0.09];
 x.TMCFreq = [250   500  1000  2000  4000  8000];
-x.TMC = [	36.8	38.9	36	37.8	48.8	45.5	 
-44	43.6	38.3	42	52.1	48	 
-49.5	51.3	49.8	50.2	62.7	53.3	 
-54	54.4	52	73.1	68.1	51.8	 
-61	71.4	56.9	73.3	72.4	57.3	 
+x.TMC = [
+37.8	32.6	28.6	40.8	41.7	38	 
+41.5	36.3	40.1	44.5	44.6	37.5	 
+48.7	42.6	43	49.4	49.7	41.2	 
+55.4	45.6	52.2	52.9	52.3	43.2	 
+62.7	55.5	62.5	52.3	58	45.9	 
 ];
 x.TMC = x.TMC';
+
 x.MaskerRatio = [0.5      0.7      0.9        1      1.1      1.3      1.6];
 x.IFMCFreq = [250   500  1000  2000  4000  8000];
-x.IFMCs = [	60.7	75.2	79.7	75.3	75.5	84.2	 
-49.2	64.1	69.6	75.9	82.2	88.1	 
-37.5	41.7	35.8	50.5	54.2	49.5	 
-37.6	38.4	37.9	37.1	53.2	45.2	 
-39.1	36.2	37.1	47.4	42.9	47.5	 
-40	39.4	55.5	84.7	61.7	87.6	 
-53.1	82.7	83.6	89.9	98.5	101	 
+x.IFMCs = [
+44.3	58.4	66.1	73.5	81.9	82.9	 
+39.6	46.8	49.9	54	66.5	67.6	 
+36.8	36	35.1	42.8	49.9	38	 
+37	32.9	33.1	40.8	41.7	37.3	 
+34.6	30.6	30.6	41.9	42	38.3	 
+35.3	31.8	42.7	52.4	61	56.7	 
+40.5	45.3	61	75.5	85.3	90.5	 
 ];
 x.IFMCs = x.IFMCs';
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/multithreshold 1.46/MTprofile10_0hr19_Aug_2011.m	Fri Aug 19 16:07:07 2011 +0100
@@ -0,0 +1,31 @@
+function x = MTprofile10_0hr19_Aug_2011
+%created: 10_0hr19_Aug_2011
+
+x.BFs = [250   500  1000  2000  4000  8000];
+
+x.LongTone = [16.1      13.1      9.24      9.67      15.4      20.6];
+x.ShortTone = [18.9        16      12.3      12.7      20.8      26.6];
+
+x.Gaps = [0.01      0.03      0.05      0.07      0.09];
+x.TMCFreq = [250   500  1000  2000  4000  8000];
+x.TMC = [
+36.2	31.2	35.2	54.4	38.2	37	 
+69.5	42.5	52.6	64.8	62.9	39.8	 
+77.8	65.2	52.2	66.7	72.3	42.2	 
+80.7	71.5	70.7	71.2	80.6	44.4	 
+89.6	82.3	77.8	74.9	78.2	68.2	 
+];
+x.TMC = x.TMC';
+
+x.MaskerRatio = [0.5      0.7      0.9        1      1.1      1.3      1.6];
+x.IFMCFreq = [250   500  1000  2000  4000  8000];
+x.IFMCs = [
+68.2	79.3	77.9	75.2	77.4	83.2	 
+58	66.9	69.2	74.7	83.1	88	 
+52.4	38.3	39.6	53.2	74.6	41.5	 
+41.7	32.1	31.3	46.9	40.4	39.4	 
+35	29.2	23.9	56.3	40.4	37.2	 
+43.3	37.7	70.7	80.1	91.7	89	 
+78.7	80.7	85.5	91	99.7	102	 
+];
+x.IFMCs = x.IFMCs';
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/multithreshold 1.46/MTprofile17_14hr18_Aug_2011.m	Fri Aug 19 16:07:07 2011 +0100
@@ -0,0 +1,31 @@
+function x = MTprofile17_14hr18_Aug_2011
+%created: 17_14hr18_Aug_2011
+
+x.BFs = [1000];
+
+x.LongTone = [9.11];
+x.ShortTone = [11.6];
+
+x.Gaps = [0.01      0.03      0.05      0.07      0.09];
+x.TMCFreq = [1000];
+x.TMC = [
+34.5	 
+37.9	 
+42	 
+52	 
+55.7	 
+];
+x.TMC = x.TMC';
+
+x.MaskerRatio = [0.5      0.7      0.9        1      1.1      1.3      1.6];
+x.IFMCFreq = [1000];
+x.IFMCs = [
+64.9	 
+48.4	 
+33.9	 
+32.7	 
+28.9	 
+38.6	 
+60.9	 
+];
+x.IFMCs = x.IFMCs';
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/multithreshold 1.46/MTprofile17_32hr18_Aug_2011.m	Fri Aug 19 16:07:07 2011 +0100
@@ -0,0 +1,8 @@
+function x = MTprofile17_32hr18_Aug_2011
+%created: 17_32hr18_Aug_2011
+
+x.BFs = [1000];
+
+x.LongTone = [9.11];
+x.ShortTone = [11.6];
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/multithreshold 1.46/MTprofile17_33hr18_Aug_2011.m	Fri Aug 19 16:07:07 2011 +0100
@@ -0,0 +1,21 @@
+function x = MTprofile17_33hr18_Aug_2011
+%created: 17_33hr18_Aug_2011
+
+x.BFs = [1000];
+
+x.LongTone = [9.11];
+x.ShortTone = [11.6];
+
+x.Gaps = [0.01      0.03      0.05      0.07      0.09];
+x.TMCFreq = [1000];
+x.TMC = [
+34.5	37.9	42	52	55.7	 
+];
+x.TMC = x.TMC';
+
+x.MaskerRatio = [0.5      0.7      0.9        1      1.1      1.3      1.6];
+x.IFMCFreq = [1000];
+x.IFMCs = [
+64.9	48.4	33.9	32.7	28.9	38.6	60.9	 
+];
+x.IFMCs = x.IFMCs';
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/multithreshold 1.46/MTprofile17_44hr18_Aug_2011.m	Fri Aug 19 16:07:07 2011 +0100
@@ -0,0 +1,31 @@
+function x = MTprofile17_44hr18_Aug_2011
+%created: 17_44hr18_Aug_2011
+
+x.BFs = [250   500  1000  2000  4000  8000];
+
+x.LongTone = [16.1      12.7      8.93      8.95        14      20.9];
+x.ShortTone = [18.8      15.5      11.9      13.6      20.2      25.3];
+
+x.Gaps = [0.01      0.03      0.05      0.07      0.09];
+x.TMCFreq = [250   500  1000  2000  4000  8000];
+x.TMC = [
+51.3	49.7	43.7	56.1	45.7	35.3	 
+59.2	52.9	57.9	63.8	58	39.2	 
+63.3	65.6	58.6	67.5	63.9	45.4	 
+77	72	70.5	68.8	69.7	39.7	 
+82.7	77.3	75.1	73.3	70.4	70.4	 
+];
+x.TMC = x.TMC';
+
+x.MaskerRatio = [0.5      0.7      0.9        1      1.1      1.3      1.6];
+x.IFMCFreq = [250   500  1000  2000  4000  8000];
+x.IFMCs = [
+61.1	74.9	78.7	76.7	78.8	84.8	 
+55.9	64.3	64.8	76.9	85.3	89.9	 
+52.4	53.1	47.8	59.3	62	37.6	 
+48.8	44	45	54.9	56.9	36.7	 
+53.7	47.4	43.2	63.6	45.9	39.7	 
+52.5	55	61.7	74.5	82.5	72.1	 
+59.2	68.7	90.1	104	102	103	 
+];
+x.IFMCs = x.IFMCs';
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/multithreshold 1.46/MTprofile7_44hr19_Aug_2011.m	Fri Aug 19 16:07:07 2011 +0100
@@ -0,0 +1,31 @@
+function x = MTprofile7_44hr19_Aug_2011
+%created: 7_44hr19_Aug_2011
+
+x.BFs = [1000];
+
+x.LongTone = [9.83];
+x.ShortTone = [12.3];
+
+x.Gaps = [0.01      0.03      0.05      0.07      0.09];
+x.TMCFreq = [1000];
+x.TMC = [
+63.5	 
+66.8	 
+72	 
+82.7	 
+90.2	 
+];
+x.TMC = x.TMC';
+
+x.MaskerRatio = [0.5      0.7      0.9        1      1.1      1.3      1.6];
+x.IFMCFreq = [1000];
+x.IFMCs = [
+81	 
+73.3	 
+63.4	 
+60.3	 
+65.5	 
+79	 
+88	 
+];
+x.IFMCs = x.IFMCs';
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/multithreshold 1.46/MTprofile8_19hr19_Aug_2011.m	Fri Aug 19 16:07:07 2011 +0100
@@ -0,0 +1,29 @@
+function x = MTprofile8_19hr19_Aug_2011
+%created: 8_19hr19_Aug_2011
+
+x.BFs = [250   500  1000  2000  4000  8000];
+
+x.LongTone = [16.5        13       9.2      9.64      15.4      21.2];
+x.ShortTone = [19.6      15.4      12.2      12.6      20.5      25.2];
+
+x.Gaps = [0.01      0.03      0.05      0.07      0.09];
+x.TMCFreq = [250   500  1000  2000  4000  8000];
+x.TMC = [
+75.2	55.2	60.6	62.6	63.2	37	 
+81	67.3	69	70.4	73.5	38.4	 
+88.3	71.1	74.2	74.6	77.7	39.4	 
+99	84.9	78.2	75.4	80.9	48.2	 
+NaN	86.9	83.4	79.8	83.8	66.2	 
+];
+x.TMC = x.TMC';
+
+x.MaskerRatio = [0.01     0.03     0.05     0.07     0.09];
+x.IFMCFreq = [250   500  1000  2000  4000  8000];
+x.IFMCs = [
+75.2	55.2	60.6	62.6	63.2	37	 
+81	67.3	69	70.4	73.5	38.4	 
+88.3	71.1	74.2	74.6	77.7	39.4	 
+99	84.9	78.2	75.4	80.9	48.2	 
+NaN	86.9	83.4	79.8	83.8	66.2	 
+];
+x.IFMCs = x.IFMCs';
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/multithreshold 1.46/MTprofile8_24hr19_Aug_2011.m	Fri Aug 19 16:07:07 2011 +0100
@@ -0,0 +1,31 @@
+function x = MTprofile8_24hr19_Aug_2011
+%created: 8_24hr19_Aug_2011
+
+x.BFs = [1000];
+
+x.LongTone = [10.7];
+x.ShortTone = [13.6];
+
+x.Gaps = [0.01      0.03      0.05      0.07      0.09];
+x.TMCFreq = [1000];
+x.TMC = [
+53.2	 
+68.6	 
+72.4	 
+75	 
+85.1	 
+];
+x.TMC = x.TMC';
+
+x.MaskerRatio = [0.5      0.7      0.9        1      1.1      1.3      1.6];
+x.IFMCFreq = [1000];
+x.IFMCs = [
+79.9	 
+73.8	 
+61	 
+52.3	 
+63.3	 
+77.8	 
+88.4	 
+];
+x.IFMCs = x.IFMCs';
--- a/multithreshold 1.46/expGUI_MT.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/multithreshold 1.46/expGUI_MT.m	Fri Aug 19 16:07:07 2011 +0100
@@ -610,19 +610,19 @@
 % ------------------------------------------------ pushbuttonRun_Callback
 function pushbuttonRun_Callback(hObject, eventdata, handles)
 global checkForPreviousGUI % holds screen positioning across repeated calls
-global experiment betweenRuns paradigmNames errormsg
+global experiment stimulusParameters betweenRuns paradigmNames errormsg
 checkForPreviousGUI.GUIposition=get(handles.figure1,'position');
 experiment.singleShot=0;
 experiment.stop=0;
 
 switch experiment.paradigm
     case 'profile'
-        %% special option for two successive and linked measurements
+        %% special option for  successive and linked measurements
+        global resultsTable
         experiment.paradigm='threshold_16ms';
         set(handles.edittargetDuration,'string', num2str(0.25))
         set(handles.editstopCriteriaBox,'string','10') % nTrials
         run (handles)
-
         if experiment.stop
             optionNo=strmatch('profile',paradigmNames);
             set(handles.popupmenuParadigm,'value',optionNo);
@@ -631,11 +631,9 @@
             return
         end
 
-        global resultsTable
         longTone=resultsTable(2:end,2:end);
-
         set(handles.edittargetDuration,'string', num2str(0.016))
-        set(handles.editstopCriteriaBox,'string','20') % nTrials
+        set(handles.editstopCriteriaBox,'string','30') % nTrials
         run (handles)
         if experiment.stop
             disp(errormsg)
@@ -646,19 +644,17 @@
             aParadigmSelection(handles)
             return
         end
-
         shortTone=resultsTable(2:end,2:end);
-
-        % use these threshold for TMC
+        
+        % use 16ms thresholds for TMC
         thresholds16ms=betweenRuns.thresholds;
         optionNo=strmatch('TMC',paradigmNames);
         set(handles.popupmenuParadigm,'value',optionNo);
         aParadigmSelection(handles)
         set(handles.edittargetLevel,'string', thresholds16ms+10);
-        set(handles.editstopCriteriaBox,'string','20')  % nTrials
+        set(handles.editstopCriteriaBox,'string','30')  % nTrials
         pause(.1)
         run (handles)
-
         if experiment.stop
             disp(errormsg)
             optionNo=strmatch('profile',paradigmNames);
@@ -668,21 +664,18 @@
             aParadigmSelection(handles)
             return
         end
-
         TMC=resultsTable(2:end,2:end);
         gaps=resultsTable(2:end,1);
         BFs=resultsTable(1, 2:end);
 
-        % use these threshold for IFMC
+        % use 16ms threshold for IFMC
         optionNo=strmatch('IFMC',paradigmNames);
         set(handles.popupmenuParadigm,'value',optionNo);
         aParadigmSelection(handles)
         set(handles.edittargetLevel,'string', thresholds16ms+10);
-        set(handles.editstopCriteriaBox,'string','10')  % nTrials
+        set(handles.editstopCriteriaBox,'string','30')  % nTrials
         pause(.1)
         run (handles)
-
-
         IFMCs=resultsTable(2:end,2:end);
         offBFs=resultsTable(2:end,1);
 
@@ -694,10 +687,11 @@
 %% save data and plot profile
 
 %         save profile longTone shortTone gaps BFs TMC offBFs IFMCs
+fileName=['MTprofile' Util_timeStamp];
 profile2mFile(longTone, shortTone, gaps, BFs, TMC, offBFs, IFMCs,...
-    'MTprofile')
-plotProfile('MTprofile', 'profile_CMA_L')
-
+    fileName)
+while ~fopen(fileName), end
+plotProfile(fileName, 'profile_CMA_L')
 %% xx
         if strcmp(errormsg,'manually stopped')
             disp(errormsg)
@@ -1068,6 +1062,7 @@
 
 % identify model parameter changes if any
 paramChanges=get(handles.editparamChanges,'string');
+if ~strcmp(paramChanges, ';'), paramChanges=[paramChanges ';']; end
 eval(paramChanges);
 
 % -------------------------------------------- aSetSampleRate
@@ -1418,6 +1413,7 @@
 relativeFrequencies=[0.25    .5   .75  1  1.25 1.5    2];
 relativeFrequencies=[ 1 ];
 AN_spikesOrProbability='probability';
+AN_spikesOrProbability='spikes';
 testBM(stimulusParameters.targetFrequency, ...
     experiment.name,relativeFrequencies, AN_spikesOrProbability, ...
     paramChanges);
@@ -1429,8 +1425,9 @@
 showPSTHs=0;
 targetFrequency=stimulusParameters.targetFrequency(1);
 BFlist=targetFrequency;
-
-testAN(targetFrequency,BFlist,-10:10:90,experiment.name, paramChanges);
+testLevels=-10:10:90;
+% testLevels=80:10:90;
+testAN(targetFrequency,BFlist,testLevels,experiment.name, paramChanges);
 
 function pushbuttonPhLk_Callback(hObject, eventdata, handles)
 global experiment
@@ -1451,7 +1448,7 @@
 aReadAndCheckParameterBoxes(handles);
 showPSTHs=1;
 testFM(stimulusParameters.targetFrequency(1),experiment.name, ...
-    'probability', paramChanges)
+    'spikes', paramChanges)
 
 function popupmenuPhase_Callback(hObject, eventdata, handles)
 global stimulusParameters
--- a/multithreshold 1.46/paradigms/paradigm_profile.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/multithreshold 1.46/paradigms/paradigm_profile.m	Fri Aug 19 16:07:07 2011 +0100
@@ -16,9 +16,7 @@
 
 stimulusParameters.rampDuration=0.004;
 
-experiment.stopCriteria2IFC=[75 3 5];
-experiment.singleIntervalMaxTrials=[20];
-
+experiment.singleIntervalMaxTrials=[50];
 
 % forced choice window interval
 stimulusParameters.AFCsilenceDuration=0.5;
--- a/multithreshold 1.46/plotProfile.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/multithreshold 1.46/plotProfile.m	Fri Aug 19 16:07:07 2011 +0100
@@ -32,7 +32,8 @@
 % TMC
 for BFno=1:length(foreground.TMCFreq)
     subplot(2,6,BFno)
-    plot(foreground.Gaps,foreground.TMC(BFno,:)-foreground.LongTone(BFno),'r','lineWidth',3), hold on
+    % SL
+% plot(foreground.Gaps,foreground.TMC(BFno,:)-foreground.LongTone(BFno),'r','lineWidth',3), hold on
     plot(foreground.Gaps,foreground.TMC(BFno,:),'b','lineWidth',3), hold on
     ylim([-10 110])
     xlim([0.01 0.1])
@@ -52,7 +53,8 @@
         if ~isempty(idx);
             
             subplot(2,6,idx)
-            plot(background.Gaps,background.TMC(BFno,:)-background.LongTone(BFno),'k:')
+            % SL
+% plot(background.Gaps,background.TMC(BFno,:)-background.LongTone(BFno),'k:')
             plot(background.Gaps,background.TMC(BFno,:),'k:')
             ylim([-10 110])
             xlim([0.01 0.1])
@@ -82,10 +84,11 @@
         xlim([100 12000])
     end
 end
-mydate=datestr(now); idx=findstr(':',mydate); mydate(idx)='_';
+title([fgName ' / ' bgName])
+% mydate=datestr(now); idx=findstr(':',mydate); mydate(idx)='_';
 
-fileName= ['savedData/' mydate ];
-
-save (fileName)
-set(gcf,'name', mydate)
-disp(fileName)
+% fileName= ['savedData/' mydate ];
+% 
+% save (fileName)
+% set(gcf,'name', mydate)
+% disp(fileName)
--- a/multithreshold 1.46/profile2mFile.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/multithreshold 1.46/profile2mFile.m	Fri Aug 19 16:07:07 2011 +0100
@@ -5,14 +5,17 @@
 fid = fopen([mFileName '.m'],'w');
 fprintf(fid, '%s\n', St);
 
+St = ['%created: ' UTIL_timeStamp];
+fprintf(fid, '%s\n\n', St);
+
 St = ['x.BFs = [' num2str(BFs) '];' ];
-fprintf(fid, '%s\n', St);
+fprintf(fid, '%s\n\n', St);
 
 St = ['x.LongTone = [' num2str(longTone',3) '];' ];
 fprintf(fid, '%s\n', St);
 
 St = ['x.ShortTone = [' num2str(shortTone',3) '];' ];
-fprintf(fid, '%s\n', St);
+fprintf(fid, '%s\n\n', St);
 
 St = ['x.Gaps = [' num2str(gaps',3) '];' ];
 fprintf(fid, '%s\n', St);
@@ -21,7 +24,7 @@
 St = ['x.TMCFreq = [' num2str(TMCFreq) '];' ];
 fprintf(fid, '%s\n', St);
 
-fprintf(fid, '%s\t', 'x.TMC = [');
+fprintf(fid, '%s\n', 'x.TMC = [');
 for i = 1:size(TMC,1),
     for j = 1:size(TMC,2),
         St = [num2str(TMC(i,j),3)];
@@ -32,7 +35,7 @@
 fprintf(fid,'%s\n','];');
 
 St = ['x.TMC = x.TMC'';' ];
-fprintf(fid, '%s\n', St);
+fprintf(fid, '%s\n\n', St);
 
 St = ['x.MaskerRatio = [' num2str(offBFs',2) '];' ];
 fprintf(fid, '%s\n', St);
@@ -41,7 +44,7 @@
 St = ['x.IFMCFreq = [' num2str(IFMCFreq) '];' ];
 fprintf(fid, '%s\n', St);
 
-fprintf(fid, '%s\t', 'x.IFMCs = [');
+fprintf(fid, '%s\n', 'x.IFMCs = [');
 for i = 1:size(IFMCs,1),
     for j = 1:size(IFMCs,2),
         St = [num2str(IFMCs(i,j),3)];
Binary file multithreshold 1.46/savedData/11-Aug-2011 10_53_52.mat has changed
Binary file multithreshold 1.46/savedData/11-Aug-2011 10_55_36.mat has changed
Binary file multithreshold 1.46/savedData/11-Aug-2011 15_51_34.mat has changed
Binary file multithreshold 1.46/savedData/11-Aug-2011 15_54_08.mat has changed
Binary file multithreshold 1.46/savedData/11-Aug-2011 15_57_29.mat has changed
Binary file multithreshold 1.46/savedData/11-Aug-2011 16_14_22.mat has changed
Binary file multithreshold 1.46/savedData/11-Aug-2011 16_26_04.mat has changed
Binary file multithreshold 1.46/savedData/11-Aug-2011 16_47_28.mat has changed
Binary file multithreshold 1.46/savedData/11-Aug-2011 16_48_53.mat has changed
Binary file multithreshold 1.46/savedData/11-Aug-2011 16_49_46.mat has changed
Binary file multithreshold 1.46/savedData/11-Aug-2011 16_50_10.mat has changed
Binary file multithreshold 1.46/savedData/11-Aug-2011 16_50_37.mat has changed
Binary file multithreshold 1.46/savedData/12-Aug-2011 05_33_35.mat has changed
Binary file multithreshold 1.46/savedData/12-Aug-2011 05_41_36.mat has changed
Binary file multithreshold 1.46/savedData/12-Aug-2011 05_48_34.mat has changed
Binary file multithreshold 1.46/savedData/12-Aug-2011 05_56_54.mat has changed
Binary file multithreshold 1.46/savedData/12-Aug-2011 06_40_17.mat has changed
Binary file multithreshold 1.46/savedData/12-Aug-2011 06_56_59.mat has changed
Binary file multithreshold 1.46/savedData/12-Aug-2011 07_20_14.mat has changed
Binary file multithreshold 1.46/savedData/12-Aug-2011 07_29_20.mat has changed
Binary file multithreshold 1.46/savedData/12-Aug-2011 07_43_17.mat has changed
Binary file multithreshold 1.46/savedData/12-Aug-2011 07_49_00.mat has changed
Binary file multithreshold 1.46/savedData/mostRecentResults.mat has changed
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/multithreshold 1.46/temp.m	Fri Aug 19 16:07:07 2011 +0100
@@ -0,0 +1,14 @@
+foreground = MTprofile17_14hr18_Aug_2011
+
+fileName=['MTprofile' Util_timeStamp];
+longTone=foreground.LongTone;
+shortTone=foreground.ShortTone;
+gaps=foreground.Gaps;
+BFs=foreground.BFs;
+TMC=foreground.TMC;
+offBFs=foreground.MaskerRatio;
+IFMCs=foreground.IFMCs;
+profile2mFile(longTone, shortTone, gaps', BFs, TMC', offBFs', IFMCs',...
+    fileName)
+pause(1)
+plotProfile(fileName, 'profile_CMA_L')
--- a/parameterStore/MAPparamsNormal.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/parameterStore/MAPparamsNormal.m	Fri Aug 19 16:07:07 2011 +0100
@@ -57,18 +57,18 @@
 % Acoustic reflex: maximum attenuation should be around 25 dB Price (1966)
 % i.e. a minimum ratio of 0.056.
 % 'spikes' model: AR based on brainstem spiking activity (LSR)
-OMEParams.rateToAttenuationFactor=0.006;   % * N(all ICspikes)
+OMEParams.rateToAttenuationFactor=0.01;   % * N(all ICspikes)
 % OMEParams.rateToAttenuationFactor=0;   % * N(all ICspikes)
 
 % 'probability model': Ar based on AN firing probabilities (LSR)
-OMEParams.rateToAttenuationFactorProb=0.01;% * N(all ANrates)
+OMEParams.rateToAttenuationFactorProb=0.006;% * N(all ANrates)
 % OMEParams.rateToAttenuationFactorProb=0;% * N(all ANrates)
 
 % asymptote should be around 100-200 ms
 OMEParams.ARtau=.05; % AR smoothing function
 % delay must be longer than the segment length
 OMEParams.ARdelay=efferentDelay;  %Moss gives 8.5 ms latency
-OMEParams.ARrateThreshold=0;
+OMEParams.ARrateThreshold=40;
 
 %%  #3 DRNL
 DRNLParams=[];  % clear the structure first
@@ -76,16 +76,17 @@
 
 % DRNL nonlinear path
 DRNLParams.a=5e4;     % DRNL.a=0 means no OHCs (no nonlinear path)
-DRNLParams.a=2e4;     % DRNL.a=0 means no OHCs (no nonlinear path)
 
 DRNLParams.b=8e-6;    % *compression threshold raised compression
 % DRNLParams.b=1;    % b=1 means no compression
+DRNLParams.cTh= 5e-8; % compression threshold in meters.
 
-DRNLParams.c=0.2;     % compression exponent
+DRNLParams.c=9e-2;     % compression exponent
 % nonlinear filters
 DRNLParams.nonlinCFs=BFlist;
 DRNLParams.nonlinOrder=	3;           % order of nonlinear gammatone filters
 p=0.2895;   q=170;  % human  (% p=0.14;   q=366;  % cat)
+p=0.2895;   q=250;  % human  (% p=0.14;   q=366;  % cat)
 DRNLParams.nlBWs=  p * BFlist + q;
 DRNLParams.p=p;   DRNLParams.q=q;   % save p and q for printing only
 
@@ -102,12 +103,12 @@
 DRNLParams.MOCdelay = efferentDelay;            % must be < segment length!
 
 % 'spikes' model: MOC based on brainstem spiking activity (HSR)
-DRNLParams.rateToAttenuationFactor = .01;  % strength of MOC
+DRNLParams.rateToAttenuationFactor = .009;  % strength of MOC
 %      DRNLParams.rateToAttenuationFactor = 0;  % strength of MOC
 % 'probability' model: MOC based on AN spiking activity (HSR)
-DRNLParams.rateToAttenuationFactorProb = .0055;  % strength of MOC
+DRNLParams.rateToAttenuationFactorProb = 0.0045;  % strength of MOC
 % DRNLParams.rateToAttenuationFactorProb = .0;  % strength of MOC
-DRNLParams.MOCrateThresholdProb =70;                % spikes/s probability only
+DRNLParams.MOCrateThresholdProb =50;                % spikes/s probability only
 
 DRNLParams.MOCtau =.1;                         % smoothing for MOC
 
@@ -148,8 +149,9 @@
 % reminder: changing z has a strong effect on HF thresholds (like Et)
 IHCpreSynapseParams.z=	    2e42;   % scalar Ca -> vesicle release rate
 
-LSRtauCa=35e-6;            HSRtauCa=85e-6;            % seconds
+LSRtauCa=35e-6;            HSRtauCa=80e-6;            % seconds
 % LSRtauCa=35e-6;            HSRtauCa=70e-6;            % seconds
+IHCpreSynapseParams.tauCa= [15e-6 80e-6]; %LSR and HSR fiber
 IHCpreSynapseParams.tauCa= [LSRtauCa HSRtauCa]; %LSR and HSR fiber
 
 %%  #6 AN_IHCsynapse
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/profiles/MTprofileMOCis0075.m	Fri Aug 19 16:07:07 2011 +0100
@@ -0,0 +1,24 @@
+function x = MTprofile
+x.BFs = [250   500  1000  2000  4000  8000];
+x.LongTone = [22.1      19.9        15      15.1      21.2      24.8];
+x.ShortTone = [25.3      21.5      17.9      18.1      24.2      29.4];
+x.Gaps = [0.01      0.03      0.05      0.07      0.09];
+x.TMCFreq = [250   500  1000  2000  4000  8000];
+x.TMC = [	40.5	43.4	41.2	47.6	42.4	42.1	 
+46.9	52.7	60.1	58.6	52.9	42.6	 
+54.3	55.5	65.6	69.3	46.8	50.1	 
+66	70.8	78.8	70.6	69.8	43.7	 
+71	91.4	80.1	75	73.6	57.8	 
+];
+x.TMC = x.TMC';
+x.MaskerRatio = [0.5      0.7      0.9        1      1.1      1.3      1.6];
+x.IFMCFreq = [250   500  1000  2000  4000  8000];
+x.IFMCs = [	61.5	77.8	82.3	75.7	75.5	79.5	 
+55.9	62.6	72.4	77.9	81	86.7	 
+39.3	46.7	50	38.9	61.5	42	 
+41.9	42.8	53	43.2	56.9	41.4	 
+41.9	40.9	42.5	44.2	53.6	41.8	 
+44.2	55.6	65	73.4	71.8	64.6	 
+85	83.2	86.9	89.1	99.6	102	 
+];
+x.IFMCs = x.IFMCs';
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/profiles/MTprofileMOCis01.m	Fri Aug 19 16:07:07 2011 +0100
@@ -0,0 +1,24 @@
+function x = MTprofile
+x.BFs = [250   500  1000  2000  4000  8000];
+x.LongTone = [22.2      18.5      13.6      14.1      19.9      25.5];
+x.ShortTone = [25.7      22.4      18.1      18.6        24      30.3];
+x.Gaps = [0.01      0.03      0.05      0.07      0.09];
+x.TMCFreq = [250   500  1000  2000  4000  8000];
+x.TMC = [	40.1	35.7	31.3	46.8	37.1	41.8	 
+41.2	41	38.7	53.3	53.5	43.2	 
+50.6	50.2	38.3	46.2	49.9	43.9	 
+55.7	54.3	53.7	62	44.7	48.6	 
+74.2	60.7	70.3	65.3	65.1	52.8	 
+];
+x.TMC = x.TMC';
+x.MaskerRatio = [0.5      0.7      0.9        1      1.1      1.3      1.6];
+x.IFMCFreq = [250   500  1000  2000  4000  8000];
+x.IFMCs = [	69.7	78.8	80	74.6	74.7	79.6	 
+48.5	64.5	70.1	76.3	80.2	88.2	 
+39.1	42.7	36.9	34.1	56.6	41.9	 
+42.1	40.3	38.2	47.4	37.9	40.1	 
+35.3	32.4	35.8	41.6	52.7	46.9	 
+39.7	39.2	49	85.4	86	61.3	 
+48.2	82.8	84.5	88.2	95.8	100	 
+];
+x.IFMCs = x.IFMCs';
--- a/testPrograms/testAN.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/testPrograms/testAN.m	Fri Aug 19 16:07:07 2011 +0100
@@ -6,62 +6,25 @@
 
 global IHC_VResp_VivoParams  IHC_cilia_RPParams IHCpreSynapseParams
 global AN_IHCsynapseParams
-
 global ANoutput ANdt CNoutput ICoutput ICmembraneOutput ANtauCas
 global ARattenuation MOCattenuation
 
 dbstop if error
 restorePath=path;
-
 addpath (['..' filesep 'MAP'], ['..' filesep 'utilities'], ...
     ['..' filesep 'parameterStore'],  ['..' filesep 'wavFileStore'],...
     ['..' filesep 'testPrograms'])
 
-if nargin<5
-    paramChanges=[];
-end
-
-if nargin<4
-    paramsName='Normal';
-end
-
+if nargin<5, paramChanges=[]; end
+if nargin<4, paramsName='Normal'; end
 if nargin<3
     levels=-10:10:80;
-    % levels=80:10:90;
 end
-
 nLevels=length(levels);
 
-toneDuration=.2;
-rampDuration=0.002;
-silenceDuration=.02;
+toneDuration=.2;   rampDuration=0.002;   silenceDuration=.02;
 localPSTHbinwidth=0.001;
 
-% Use only the first frequency in the GUI targetFrequency box to defineBF
-% targetFrequency=stimulusParameters.targetFrequency(1);
-% BFlist=targetFrequency;
-
-AN_HSRonset=zeros(nLevels,1);
-AN_HSRsaturated=zeros(nLevels,1);
-AN_LSRonset=zeros(nLevels,1);
-AN_LSRsaturated=zeros(nLevels,1);
-CNLSRrate=zeros(nLevels,1);
-CNHSRsaturated=zeros(nLevels,1);
-ICHSRsaturated=zeros(nLevels,1);
-ICLSRsaturated=zeros(nLevels,1);
-vectorStrength=zeros(nLevels,1);
-
-AR=zeros(nLevels,1);
-MOC=zeros(nLevels,1);
-
-% ANoutput=zeros(200,200);
-
-figure(15), clf
-set(gcf,'position',[980   356   401   321])
-figure(5), clf
-set(gcf,'position', [980 34 400 295])
-drawnow
-
 %% guarantee that the sample rate is at least 10 times the frequency
 sampleRate=50000;
 while sampleRate< 10* targetFrequency
@@ -77,16 +40,36 @@
 ANdt=period/round(period/ANdt);
 dt=ANdt/ANspeedUpFactor;
 sampleRate=1/dt;
-epochsPerPeriod=sampleRate*period;
+
+%% pre-allocate storage
+AN_HSRonset=zeros(nLevels,1);
+AN_HSRsaturated=zeros(nLevels,1);
+AN_LSRonset=zeros(nLevels,1);
+AN_LSRsaturated=zeros(nLevels,1);
+CNLSRrate=zeros(nLevels,1);
+CNHSRsaturated=zeros(nLevels,1);
+ICHSRsaturated=zeros(nLevels,1);
+ICLSRsaturated=zeros(nLevels,1);
+vectorStrength=zeros(nLevels,1);
+
+AR=zeros(nLevels,1);
+MOC=zeros(nLevels,1);
+
+figure(15), clf
+set(gcf,'position',[980   356   401   321])
+figure(5), clf
+set(gcf,'position', [980 34 400 295])
+drawnow
+
 
 %% main computational loop (vary level)
 levelNo=0;
 for leveldB=levels
     levelNo=levelNo+1;
-    
+    amp=28e-6*10^(leveldB/20);
     fprintf('%4.0f\t', leveldB)
-    amp=28e-6*10^(leveldB/20);
-    
+
+    %% generate tone and silences
     time=dt:dt:toneDuration;
     rampTime=dt:dt:rampDuration;
     ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
@@ -95,24 +78,20 @@
     
     silence=zeros(1,round(silenceDuration/dt));
     
-    % create signal (leveldB/ targetFrequency)
     inputSignal=amp*sin(2*pi*targetFrequency'*time);
     inputSignal= ramp.*inputSignal;
     inputSignal=[silence inputSignal];
     
     %% run the model
-    AN_spikesOrProbability='spikes';
-    showPlotsAndDetails=0;
-    
-    
+    AN_spikesOrProbability='spikes'    
     MAP1_14(inputSignal, 1/dt, BFlist, ...
         paramsName, AN_spikesOrProbability, paramChanges);
     
     nTaus=length(ANtauCas);
     
+    %% Auditory nerve evaluate and display (Fig. 5)
     %LSR (same as HSR if no LSR fibers present)
     [nANFibers nTimePoints]=size(ANoutput);
-    
     numLSRfibers=nANFibers/nTaus;
     numHSRfibers=numLSRfibers;
     
@@ -136,6 +115,7 @@
     hold on,  bar(PSTHtime,PSTHLSR,'r')
     ylim([0 1000])
     xlim([0 length(PSTH)*localPSTHbinwidth])
+    grid on
     set(gcf,'name',['PSTH: ' num2str(BFlist), ' Hz: ' num2str(leveldB) ' dB']);
     
     % AN - CV
@@ -153,9 +133,9 @@
     title(['AN HSR: CV=' num2str(cvANHSR(3),'%5.2f') ...
         'VS=' num2str(VS,'%5.2f')])
     
-    % CN - first-order neurons
+    %% CN - first-order neurons
     
-    % CN LSR
+    % CN driven by LSR fibers
     [nCNneurons c]=size(CNoutput);
     nLSRneurons=round(nCNneurons/nTaus);
     CNLSRspikes=CNoutput(1:nLSRneurons,:);
@@ -179,7 +159,7 @@
     cvMMHSR= UTIL_CV(MacGregorMultiHSRspikes, ANdt);
     title(['CN    CV= ' num2str(cvMMHSR(3),'%5.2f')])
     
-    % IC LSR
+    %% IC LSR
     [nICneurons c]=size(ICoutput);
     nLSRneurons=round(nICneurons/nTaus);
     ICLSRspikes=ICoutput(1:nLSRneurons,:);
@@ -188,18 +168,16 @@
     
     %IC HSR
     MacGregorMultiHSRspikes=...
-        ICoutput(end-nLSRneurons:end,:);
+        ICoutput(end-nLSRneurons+1:end,:);
     PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, ANdt, localPSTHbinwidth);
-    PSTH=sum(PSTH)/nLSRneurons;
-    PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only)
-    
-    ICHSRsaturated(levelNo)=mean(PSTH(length(PSTH)/2:end));
-    
+    ICHSRsaturated(levelNo)= (sum(PSTH)/nLSRneurons)/toneDuration;
+
     AR(levelNo)=min(ARattenuation);
     MOC(levelNo)=min(MOCattenuation(length(MOCattenuation)/2:end));
     
     time=dt:dt:dt*size(ICmembraneOutput,2);
     figure(5), subplot(2,2,4)
+    % plot HSR (last of two)
     plot(time,ICmembraneOutput(2, 1:end),'k')
     ylim([-0.07 0])
     xlim([0 max(time)])
@@ -207,14 +185,18 @@
     drawnow
     
     figure(5), subplot(2,2,1)
-    plot(20*log10(MOC), 'k'),
+    plot(20*log10(MOC), 'k'), hold on
+    plot(20*log10(AR), 'r'), hold off
     title(' MOC'), ylabel('dB attenuation')
     ylim([-30 0])
     
     
 end % level
+
+%% plot with levels on x-axis
 figure(5), subplot(2,2,1)
-plot(levels,20*log10(MOC), 'k'),
+plot(levels,20*log10(MOC), 'k'),hold on
+plot(levels,20*log10(AR), 'r'), hold off
 title(' MOC'), ylabel('dB attenuation')
 ylim([-30 0])
 xlim([0 max(levels)])
@@ -230,9 +212,7 @@
 fprintf('%6.2f\t', levels)
 fprintf('\n')
 
-
-% ---------------------------------------------------- display parameters
-
+%% ---------------------- display summary results (Fig 15)
 figure(15), clf
 nRows=2; nCols=2;
 
@@ -244,7 +224,7 @@
 ttl=['tauCa= ' num2str(IHCpreSynapseParams.tauCa)];
 title( ttl)
 xlabel('level dB SPL'), ylabel('peak rate (sp/s)'), grid on
-text(0, 800, 'AN onset', 'fontsize', 16)
+text(0, 800, 'AN onset', 'fontsize', 14)
 
 % AN rate - level ADAPTED function
 subplot(nRows,nCols,2)
@@ -259,7 +239,7 @@
     '  sat=' num2str(mean(AN_HSRsaturated(end,1)),'%4.0f')];
 title( ttl)
 xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')
-text(0, 340, 'AN adapted', 'fontsize', 16), grid on
+text(0, 340, 'AN adapted', 'fontsize', 14), grid on
 
 % CN rate - level ADAPTED function
 subplot(nRows,nCols,3)
@@ -274,7 +254,7 @@
     num2str(mean(CNHSRsaturated(end,1)),'%4.0f')];
 title( ttl)
 xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')
-text(0, 350, 'CN', 'fontsize', 16), grid on
+text(0, 350, 'CN', 'fontsize', 14), grid on
 
 % IC rate - level ADAPTED function
 subplot(nRows,nCols,4)
@@ -284,17 +264,13 @@
 set(gca,'ytick',0:50:300)
 xlim([min(levels) max(levels)])
 set(gca,'xtick',[levels(1):20:levels(end)]), grid on
-
-
 ttl=['spont=' num2str(mean(ICHSRsaturated(1,:)),'%4.0f') ...
     '  sat=' num2str(mean(ICHSRsaturated(end,1)),'%4.0f')];
 title( ttl)
 xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')
-text(0, 350, 'IC', 'fontsize', 16)
+text(0, 350, 'IC', 'fontsize', 14)
 set(gcf,'name',' AN CN IC rate/level')
 
-peakVectorStrength=max(vectorStrength);
-
 fprintf('\n')
 disp('levels vectorStrength')
 fprintf('%3.0f \t %6.4f \n', [levels; vectorStrength'])
--- a/testPrograms/testANprob.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/testPrograms/testANprob.m	Fri Aug 19 16:07:07 2011 +0100
@@ -1,6 +1,8 @@
-function vectorStrength=testANprob(targetFrequency,BFlist, levels, ...
+function testANprob(targetFrequency,BFlist, levels, ...
     paramsName, paramChanges)
 
+% testANprob(1000,1000, -10:10:80, 'Normal')
+
 global IHC_VResp_VivoParams  IHC_cilia_RPParams IHCpreSynapseParams
 global AN_IHCsynapseParams
 global ANprobRateOutput dt ANtauCas
@@ -45,7 +47,7 @@
 MOC=zeros(nLevels,1);
 
 figure(15), clf
-set(gcf,'position',[607    17   368   321])
+set(gcf,'position',[980   356   401   321])
 drawnow
 
 %% guarantee that the sample rate is at least 10 times the frequency
@@ -122,15 +124,17 @@
 
 
     figure(15), subplot(2,2,3)
-    plot(20*log10(MOC), 'k'),
-    title(' MOC'), ylabel('dB attenuation')
+    plot(20*log10(MOC), 'k'), hold on
+    plot(20*log10(AR), 'r'),  hold off
+    title(' MOC/AR'), ylabel('dB attenuation')
     ylim([-30 0])
 
+end % level
 
-end % level
 figure(15), subplot(2,2,3)
-plot(levels,20*log10(MOC), 'k'),
-title(' MOC'), ylabel('dB attenuation')
+plot(levels,20*log10(MOC), 'k'), hold on
+plot(levels,20*log10(AR), 'r'),  hold off
+title(' MOC/AR'), ylabel('dB attenuation')
 ylim([-30 0])
 xlim([0 max(levels)])
 
@@ -159,13 +163,13 @@
 ttl=['tauCa= ' num2str(IHCpreSynapseParams.tauCa)];
 title( ttl)
 xlabel('level dB SPL'), ylabel('peak rate (sp/s)'), grid on
-text(0, 800, 'AN onset', 'fontsize', 16)
+text(0, 800, 'AN onset', 'fontsize', 14)
 
 % AN rate - level ADAPTED function
 subplot(nRows,nCols,2)
 plot(levels,AN_LSRsaturated, 'ro'), hold on
 plot(levels,AN_HSRsaturated, 'ko'), hold off
-ylim([0 400])
+ylim([0 340])
 set(gca,'ytick',0:50:300)
 xlim([min(levels) max(levels)])
 set(gca,'xtick',[levels(1):20:levels(end)])
@@ -174,7 +178,7 @@
     '  sat=' num2str(mean(AN_HSRsaturated(end,1)),'%4.0f')];
 title( ttl)
 xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')
-text(0, 340, 'AN adapted', 'fontsize', 16), grid on
+text(0, 340, 'AN adapted', 'fontsize', 14), grid on
 
 allData=[ levels'  AN_HSRonset AN_HSRsaturated...
     AN_LSRonset AN_LSRsaturated ];
@@ -186,8 +190,6 @@
 UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams')
 UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams')
 
-fprintf('\n')
-disp('levels vectorStrength')
 
 allData=[ levels'  AN_HSRonset AN_HSRsaturated...
     AN_LSRonset AN_LSRsaturated ];
--- a/testPrograms/testBM.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/testPrograms/testBM.m	Fri Aug 19 16:07:07 2011 +0100
@@ -9,6 +9,7 @@
 % Tuning curves are generated based on a range of frequencies reletove to
 % the BF of the location.
 %
+% testBM (1000, 'Normal', 1, 'probability', [])
 
 global    DRNLParams
 
@@ -17,7 +18,7 @@
 end
 
 if nargin<4
-    AN_spikesOrProbability='probability';
+    AN_spikesOrProbability='spikes';
 end
 
 savePath=path;
@@ -165,5 +166,10 @@
 
     UTIL_printTabTable([levels' finalSummary], ...
         num2str([0 stimulusFrequencies]','%5.0f'), '%5.0f')
+    diff(finalSummary)
+    if ~isempty(paramChanges)
+        disp(paramChanges)
+    end
+    
 
 path(savePath);
--- a/testPrograms/testFM.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/testPrograms/testFM.m	Fri Aug 19 16:07:07 2011 +0100
@@ -144,7 +144,8 @@
 
         % **********************************  run MAP model
 %         showPlotsAndDetails=0;
-
+nChanges=length(paramChanges);
+paramChanges{nChanges+1}='AN_IHCsynapseParams.numFibers=	500;'; 
         MAP1_14(inputSignal, 1/dt, targetFrequency, ...
             paramsName, AN_spikesOrProbability, paramChanges);
 
@@ -159,7 +160,7 @@
             ANresponse=ANoutput(end-nLSRfibers:end,:);
         end
 
-        ANresponse=sum(ANresponse)/nLSRfibers;
+        ANresponse=sum(ANresponse);
 %         ANresponseTimes=ANdt:ANdt:length(ANresponse)*ANdt;
 %         figure(99), plot(ANresponseTimes,ANresponse)
 
@@ -189,10 +190,11 @@
             subplot(nLevels,nDurations,PSTHplotCount)
             PSTHtime=PSTHbinWidth:PSTHbinWidth:...
                 PSTHbinWidth*length(PSTH);
+            if strcmp(AN_spikesOrProbability, 'spikes')
+            bar(PSTHtime, PSTH/PSTHbinWidth/nFibers)
+%                 ylim([0 500])
+            else
             bar(PSTHtime, PSTH)
-            if strcmp(AN_spikesOrProbability, 'spikes')
-                ylim([0 500])
-            else
                 ylim([0 500])
             end
 %             xlim([0 longestSignalDuration])
--- a/testPrograms/testPhysiology.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/testPrograms/testPhysiology.m	Fri Aug 19 16:07:07 2011 +0100
@@ -1,6 +1,6 @@
 function testPhysiology(BF,paramsName, paramChanges)
 % e.g.
-% testPhysiology(2000,'Normal', [])
+% testPhysiology(1000,'Normal', [])
 
 restorePath=path;
 addpath (['..' filesep 'MAP'])
--- a/testPrograms/testPhysiologyProb.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/testPrograms/testPhysiologyProb.m	Fri Aug 19 16:07:07 2011 +0100
@@ -1,6 +1,6 @@
 function testPhysiologyProb(BF,paramsName, paramChanges)
 % e.g.
-% testPhysiologyProb(2000,'Normal', [])
+% testPhysiologyProb(1000,'Normal', [])
 
 restorePath=path;
 addpath (['..' filesep 'MAP'])
--- a/testPrograms/test_MAP1_14.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/testPrograms/test_MAP1_14.m	Fri Aug 19 16:07:07 2011 +0100
@@ -45,7 +45,6 @@
 AN_spikesOrProbability='spikes';
 
 % or
-% NB probabilities are not corrected for refractory effects
 % AN_spikesOrProbability='probability';
 
 
--- a/utilities/Util_timeStamp.m	Thu Aug 11 16:52:25 2011 +0100
+++ b/utilities/Util_timeStamp.m	Fri Aug 19 16:07:07 2011 +0100
@@ -1,4 +1,4 @@
-function timeStamp=Util_timeStamp
+function timeStamp=UTIL_timeStamp
 
 % returns a string showing time now and current date
 % the string should be suitable for creating fileNames