# HG changeset patch
# User Ray Meddis <rmeddis@essex.ac.uk>
# Date 1307985665 -3600
# Node ID 45f28c49461e7ff6cf01b53be4176102b4f47623
# Parent  c489ebada16e4d16bd7a68d38982bcbd4e4cebe8# Parent  fafe69c43108976a4409ca95e2e414d5f6a88571
removing duplicate changes

diff -r fafe69c43108 -r 45f28c49461e old files/MAPparamsNormal.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/old files/MAPparamsNormal.m	Mon Jun 13 18:21:05 2011 +0100
@@ -0,0 +1,350 @@
+function method=MAPparamsNormal ...
+    (BFlist, sampleRate, showParams)
+% MAPparams<> establishes a complete set of MAP parameters
+% Parameter file names must be of the form <MAPparams> <name>
+%
+% input arguments
+%  BFlist     (optional) specifies the desired list of channel BFs
+%    otherwise defaults set below
+%  sampleRate (optional), default is 50000.
+%  showParams (optional) =1 prints out the complete set of parameters
+% output argument
+%  method passes a miscelleny of values
+
+global inputStimulusParams OMEParams DRNLParams
+global IHC_VResp_VivoParams IHCpreSynapseParams  AN_IHCsynapseParams
+global MacGregorParams MacGregorMultiParams  filteredSACFParams
+global experiment % used by calls from multiThreshold only
+global IHC_cilia_RPParams
+
+currentFile=mfilename;                      % i.e. the name of this mfile
+method.parameterSource=currentFile(10:end); % for the record
+
+efferentDelay=0.010;
+method.segmentDuration=efferentDelay;
+
+if nargin<3, showParams=0; end
+if nargin<2, sampleRate=50000; end
+if nargin<1 || BFlist(1)<0 % if BFlist= -1, set BFlist to default
+    lowestBF=250; 	highestBF= 8000; 	numChannels=21;
+    % 21 chs (250-8k)includes BFs at 250 500 1000 2000 4000 8000
+    BFlist=round(logspace(log10(lowestBF),log10(highestBF),numChannels));
+end
+% BFlist=1000;
+
+% preserve for backward campatibility
+method.nonlinCF=BFlist; 
+method.dt=1/sampleRate; 
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% set  model parameters
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%  #1 inputStimulus
+inputStimulusParams=[];
+inputStimulusParams.sampleRate= sampleRate; 
+
+%%  #2 outerMiddleEar
+OMEParams=[];  % clear the structure first
+% outer ear resonances band pass filter  [gain lp order  hp]
+OMEParams.externalResonanceFilters=      [ 10 1 1000 4000];
+                                       
+% highpass stapes filter  
+%  Huber gives 2e-9 m at 80 dB and 1 kHz (2e-13 at 0 dB SPL)
+OMEParams.OMEstapesLPcutoff= 1000;
+OMEParams.stapesScalar=	     45e-9;
+
+% 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.004;   % * N(all ICspikes)
+%     OMEParams.rateToAttenuationFactor=0;   % * N(all ICspikes)
+
+% 'probability model': Ar based on AN firing probabilities (LSR)
+OMEParams.rateToAttenuationFactorProb=0.003;% * 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;
+
+%%  #3 DRNL
+DRNLParams=[];  % clear the structure first
+DRNLParams.BFlist=BFlist;
+
+% DRNL nonlinear path
+DRNLParams.a=3e4;     % nonlinear path gain (below compression threshold)
+<<<<<<< HEAD
+DRNLParams.a=3e2;     % DRNL.a=0 means no OHCs (no nonlinear path)
+=======
+DRNLParams.a=5e2;     % DRNL.a=0 means no OHCs (no nonlinear path)
+>>>>>>> 77bdf847c4be95767d3ea519d77878f9414a2686
+
+DRNLParams.b=8e-6;    % *compression threshold raised compression
+% DRNLParams.b=1;    % b=1 means no compression
+
+DRNLParams.c=0.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)
+DRNLParams.nlBWs=  p * BFlist + q;
+DRNLParams.p=p;   DRNLParams.q=q;   % save p and q for printing only
+
+% DRNL linear path:
+DRNLParams.g=100;     % linear path gain factor
+% linCF is not necessarily the same as nonlinCF
+minLinCF=153.13; coeffLinCF=0.7341;   % linCF>nonlinBF for BF < 1 kHz
+DRNLParams.linCFs=minLinCF+coeffLinCF*BFlist;
+DRNLParams.linOrder=	3;           % order of linear gammatone filters
+minLinBW=100; coeffLinBW=0.6531;
+DRNLParams.linBWs=minLinBW + coeffLinBW*BFlist; % bandwidths of linear  filters
+
+% DRNL MOC efferents
+DRNLParams.MOCdelay = efferentDelay;            % must be < segment length!
+% 'spikes' model: MOC based on brainstem spiking activity (HSR)
+DRNLParams.rateToAttenuationFactor = .009;  % strength of MOC
+DRNLParams.rateToAttenuationFactor = .004;  % strength of MOC
+%      DRNLParams.rateToAttenuationFactor = 0;  % strength of MOC
+
+% 'probability' model: MOC based on AN spiking activity (HSR)
+<<<<<<< HEAD
+DRNLParams.rateToAttenuationFactorProb = .007;  % strength of MOC
+DRNLParams.rateToAttenuationFactorProb = .002;  % strength of MOC
+% DRNLParams.rateToAttenuationFactorProb = .0;  % strength of MOC
+DRNLParams.MOCtau =.03;                         % smoothing for MOC
+=======
+DRNLParams.rateToAttenuationFactorProb = .004;  % strength of MOC
+% DRNLParams.rateToAttenuationFactorProb = .0;  % strength of MOC
+DRNLParams.MOCtau =.1;                         % smoothing for MOC
+>>>>>>> 77bdf847c4be95767d3ea519d77878f9414a2686
+DRNLParams.MOCrateThreshold =50;                % set to AN rate threshold
+
+
+%% #4 IHC_cilia_RPParams
+
+IHC_cilia_RPParams.tc=	0.0003;   % 0.0003 filter time simulates viscocity
+% IHC_cilia_RPParams.tc=	0.0005;   % 0.0003 filter time simulates viscocity
+IHC_cilia_RPParams.C=	0.05;      % 0.1 scalar (C_cilia ) 
+IHC_cilia_RPParams.u0=	5e-9;       
+IHC_cilia_RPParams.s0=	30e-9;
+IHC_cilia_RPParams.u1=	1e-9;
+IHC_cilia_RPParams.s1=	1e-9;
+
+IHC_cilia_RPParams.Gmax= 5e-9;    % 2.5e-9 maximum conductance (Siemens)
+IHC_cilia_RPParams.Ga=	1e-9;  % 4.3e-9 fixed apical membrane conductance
+
+%  #5 IHC_RP
+IHC_cilia_RPParams.Cab=	4e-012;         % IHC capacitance (F)
+IHC_cilia_RPParams.Cab=	1e-012;         % IHC capacitance (F)
+IHC_cilia_RPParams.Et=	0.100;          % endocochlear potential (V)
+
+IHC_cilia_RPParams.Gk=	2e-008;         % 1e-8 potassium conductance (S)
+IHC_cilia_RPParams.Ek=	-0.08;          % -0.084 K equilibrium potential
+IHC_cilia_RPParams.Rpc=	0.04;           % combined resistances
+
+
+%%  #5 IHCpreSynapse
+IHCpreSynapseParams=[];
+IHCpreSynapseParams.GmaxCa=	14e-9;% maximum calcium conductance
+IHCpreSynapseParams.GmaxCa=	12e-9;% maximum calcium conductance
+IHCpreSynapseParams.ECa=	0.066;  % calcium equilibrium potential
+IHCpreSynapseParams.beta=	400;	% determine Ca channel opening
+IHCpreSynapseParams.gamma=	100;	% determine Ca channel opening
+IHCpreSynapseParams.tauM=	0.00005;	% membrane time constant ?0.1ms
+IHCpreSynapseParams.power=	3;
+% 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=70e-6;            % seconds
+IHCpreSynapseParams.tauCa= [LSRtauCa HSRtauCa]; %LSR and HSR fiber
+
+%%  #6 AN_IHCsynapse
+% c=kym/(y(l+r)+kl)	(spontaneous rate)
+% c=(approx)  ym/l  (saturated rate)
+AN_IHCsynapseParams=[];             % clear the structure first
+AN_IHCsynapseParams.M=	12;         % maximum vesicles at synapse
+AN_IHCsynapseParams.y=	4;          % depleted vesicle replacement rate
+AN_IHCsynapseParams.y=	6;          % depleted vesicle replacement rate
+
+AN_IHCsynapseParams.x=	30;         % replenishment from re-uptake store
+AN_IHCsynapseParams.x=	60;         % replenishment from re-uptake store
+
+% reduce l to increase saturated rate
+AN_IHCsynapseParams.l=	100; % *loss rate of vesicles from the cleft
+AN_IHCsynapseParams.l=	250; % *loss rate of vesicles from the cleft
+
+AN_IHCsynapseParams.r=	500; % *reuptake rate from cleft into cell
+% AN_IHCsynapseParams.r=	300; % *reuptake rate from cleft into cell
+
+AN_IHCsynapseParams.refractory_period=	0.00075;
+% number of AN fibers at each BF (used only for spike generation)
+AN_IHCsynapseParams.numFibers=	100; 
+AN_IHCsynapseParams.TWdelay=0.004;  % ?delay before stimulus first spike
+
+AN_IHCsynapseParams.ANspeedUpFactor=5; % longer epochs for computing spikes.
+
+%%  #7 MacGregorMulti (first order brainstem neurons)
+MacGregorMultiParams=[];
+MacGregorMultiType='chopper'; % MacGregorMultiType='primary-like'; %choose
+switch MacGregorMultiType
+    case 'primary-like'
+        MacGregorMultiParams.nNeuronsPerBF=	10;   % N neurons per BF
+        MacGregorMultiParams.type = 'primary-like cell';
+        MacGregorMultiParams.fibersPerNeuron=4;   % N input fibers
+        MacGregorMultiParams.dendriteLPfreq=200;  % dendritic filter
+        MacGregorMultiParams.currentPerSpike=0.11e-6; % (A) per spike
+        MacGregorMultiParams.Cap=4.55e-9;   % cell capacitance (Siemens)
+        MacGregorMultiParams.tauM=5e-4;     % membrane time constant (s)
+        MacGregorMultiParams.Ek=-0.01;      % K+ eq. potential (V)
+        MacGregorMultiParams.dGkSpike=3.64e-5; % K+ cond.shift on spike,S
+        MacGregorMultiParams.tauGk=	0.0012; % K+ conductance tau (s)
+        MacGregorMultiParams.Th0=	0.01;   % equilibrium threshold (V)
+        MacGregorMultiParams.c=	0.01;       % threshold shift on spike, (V)
+        MacGregorMultiParams.tauTh=	0.015;  % variable threshold tau
+        MacGregorMultiParams.Er=-0.06;      % resting potential (V)
+        MacGregorMultiParams.Eb=0.06;       % spike height (V)
+
+    case 'chopper'
+        MacGregorMultiParams.nNeuronsPerBF=	10;   % N neurons per BF
+        MacGregorMultiParams.type = 'chopper cell';
+        MacGregorMultiParams.fibersPerNeuron=10;  % N input fibers
+%         MacGregorMultiParams.fibersPerNeuron=6;  % N input fibers
+
+        MacGregorMultiParams.dendriteLPfreq=50;   % dendritic filter
+        MacGregorMultiParams.currentPerSpike=35e-9; % *per spike
+        MacGregorMultiParams.currentPerSpike=30e-9; % *per spike
+        
+        MacGregorMultiParams.Cap=1.67e-8; % ??cell capacitance (Siemens)
+        MacGregorMultiParams.tauM=0.002;  % membrane time constant (s)
+        MacGregorMultiParams.Ek=-0.01;    % K+ eq. potential (V)
+        MacGregorMultiParams.dGkSpike=1.33e-4; % K+ cond.shift on spike,S
+        MacGregorMultiParams.tauGk=	0.0005;% K+ conductance tau (s)
+        MacGregorMultiParams.Th0=	0.01; % equilibrium threshold (V)
+        MacGregorMultiParams.c=	0;        % threshold shift on spike, (V)
+        MacGregorMultiParams.tauTh=	0.02; % variable threshold tau
+        MacGregorMultiParams.Er=-0.06;    % resting potential (V)
+        MacGregorMultiParams.Eb=0.06;     % spike height (V)
+        MacGregorMultiParams.PSTHbinWidth=	1e-4;
+end
+
+%%  #8 MacGregor (second-order neuron). Only one per channel
+MacGregorParams=[];                 % clear the structure first
+MacGregorParams.type = 'chopper cell';
+MacGregorParams.fibersPerNeuron=10; % N input fibers
+MacGregorParams.dendriteLPfreq=100; % dendritic filter
+MacGregorParams.currentPerSpike=120e-9;% *(A) per spike
+MacGregorParams.currentPerSpike=30e-9;% *(A) per spike
+
+MacGregorParams.Cap=16.7e-9;        % cell capacitance (Siemens)
+MacGregorParams.tauM=0.002;         % membrane time constant (s)
+MacGregorParams.Ek=-0.01;           % K+ eq. potential (V)
+MacGregorParams.dGkSpike=1.33e-4;   % K+ cond.shift on spike,S
+MacGregorParams.tauGk=	0.0005;     % K+ conductance tau (s)
+MacGregorParams.Th0=	0.01;       % equilibrium threshold (V)
+MacGregorParams.c=	0;              % threshold shift on spike, (V)
+MacGregorParams.tauTh=	0.02;       % variable threshold tau
+MacGregorParams.Er=-0.06;           % resting potential (V)
+MacGregorParams.Eb=0.06;            % spike height (V)
+MacGregorParams.debugging=0;        % (special)
+% wideband accepts input from all channels (of same fiber type)
+% use wideband to create inhibitory units
+MacGregorParams.wideband=0;         % special for wideband units
+% MacGregorParams.saveAllData=0;
+
+%%  #9 filteredSACF
+minPitch=	300; maxPitch=	3000; numPitches=60;    % specify lags
+pitches=100*log10(logspace(minPitch/100, maxPitch/100, numPitches));
+filteredSACFParams.lags=1./pitches;     % autocorrelation lags vector
+filteredSACFParams.acfTau=	.003;       % time constant of running ACF
+filteredSACFParams.lambda=	0.12;       % slower filter to smooth ACF
+filteredSACFParams.plotFilteredSACF=1;  % 0 plots unfiltered ACFs
+filteredSACFParams.plotACFs=0;          % special plot (see code)
+%  filteredSACFParams.usePressnitzer=0; % attenuates ACF at  long lags
+filteredSACFParams.lagsProcedure=  'useAllLags';
+% filteredSACFParams.lagsProcedure=  'useBernsteinLagWeights';
+% filteredSACFParams.lagsProcedure=  'omitShortLags';
+filteredSACFParams.criterionForOmittingLags=3;
+
+% checks
+if AN_IHCsynapseParams.numFibers<MacGregorMultiParams.fibersPerNeuron
+    error('MacGregorMulti: too few input fibers for input to MacG unit')
+end
+
+
+%% write all parameters to the command window
+% showParams is currently set at the top of htis function
+if showParams
+    fprintf('\n %%%%%%%%\n')
+    fprintf('\n%s\n', method.parameterSource)
+    fprintf('\n')
+    nm=UTIL_paramsList(whos);
+    for i=1:length(nm)
+        %         eval(['UTIL_showStruct(' nm{i} ', ''' nm{i} ''')'])
+        if ~strcmp(nm(i), 'method')
+            eval(['UTIL_showStructureSummary(' nm{i} ', ''' nm{i} ''', 10)'])
+        end
+    end
+end
+
+
+
+% **********************************************************************  comparison data
+% store individual data here for display on the multiThreshold GUI (if used)
+% the final value in each vector is an identifier (BF or duration))
+if isstruct(experiment)
+    switch experiment.paradigm
+        case {'IFMC','IFMC_8ms'}
+            % based on MPa
+            comparisonData=[
+                66	51	49	48	46	45	54	250;
+                60	54	46	42	39	49	65	500;
+                64	51	38	32	33	59	75	1000;
+                59	51	36	30	41	81	93	2000;
+                71	63	53	44	36	76	95	4000;
+                70	64	43	35	35	66	88	6000;
+                110	110	110	110	110	110	110	8000;
+                ];
+            if length(BFlist)==1 && ~isempty(comparisonData)
+                availableFrequencies=comparisonData(:,end)';
+                findRow= find(BFlist==availableFrequencies);
+                if ~isempty (findRow)
+                    experiment.comparisonData=comparisonData(findRow,:);
+                end
+            end
+
+        case {'TMC','TMC_8ms'}
+            % based on MPa
+            comparisonData=[
+                48	58	63	68	75	80	85	92	99	250;
+                33	39	40	49	52	61	64	77	79	500;
+                39	42	50	81	83	92	96	97	110	1000;
+                24	26	32	37	46	51	59	71	78	2000;
+                65	68	77	85	91	93	110	110	110	4000;
+                20	19	26	44	80	95	96	110	110	6000;
+                ];
+            if length(BFlist)==1 && ~isempty(comparisonData)
+                availableFrequencies=comparisonData(:,end)';
+                findRow= find(BFlist==availableFrequencies);
+                if ~isempty (findRow)
+                    experiment.comparisonData=comparisonData(findRow,:);
+                end
+            end
+
+        case { 'absThreshold',  'absThreshold_8'}
+            % MPa thresholds
+            experiment.comparisonData=[
+                32	26	16	18	22	22 0.008;
+                16	13	6	9	15	11 0.500
+                ];
+
+
+        otherwise
+            experiment.comparisonData=[];
+    end
+end
+
+
diff -r fafe69c43108 -r 45f28c49461e old files/showMAP.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/old files/showMAP.m	Mon Jun 13 18:21:05 2011 +0100
@@ -0,0 +1,336 @@
+function showMAP ()
+% defaults
+% options.showModelParameters=1;
+% options.showModelOutput=1;
+% options.printFiringRates=1;
+% options.showACF=1;
+% options.showEfferent=1;
+% options.surfProbability=0;
+% options.fileName=[];
+
+dbstop if warning
+
+global dt ANdt saveAN_spikesOrProbability savedBFlist saveMAPparamsName...
+    savedInputSignal TMoutput OMEoutput ARattenuation ...
+    DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
+    IHCoutput ANprobRateOutput ANoutput savePavailable tauCas  ...
+    CNoutput  ICoutput ICmembraneOutput ICfiberTypeRates MOCattenuation
+global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams
+global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams
+global showMapOptions
+
+
+restorePath=path;
+addpath ( ['..' filesep 'utilities'], ['..' filesep 'parameterStore'])
+
+if nargin<1
+    options.showModelParameters=1;
+    options.showModelOutput=1;
+    options.printFiringRates=1;
+    options.showACF=0;
+    options.showEfferent=1;
+    options.surfProbability=0;
+    options.fileName=[];
+end
+
+if options.showModelParameters
+    % Read parameters from MAPparams<***> file in 'parameterStore' folder
+    %  and print out all parameters
+    cmd=['MAPparams' saveMAPparamsName ...
+        '(-1, 1/dt, 1);'];
+    eval(cmd);
+end
+
+if options.printFiringRates
+    %% print summary firing rates
+    fprintf('\n\n')
+    disp('summary')
+    disp(['AR: ' num2str(min(ARattenuation))])
+    disp(['MOC: ' num2str(min(min(MOCattenuation)))])
+    nANfiberTypes=length(tauCas);
+    if strcmp(saveAN_spikesOrProbability, 'spikes')
+        nANfibers=size(ANoutput,1);
+        nHSRfibers=nANfibers/nANfiberTypes;
+        duration=size(TMoutput,2)*dt;
+        disp(['AN: ' num2str(sum(sum(ANoutput(end-nHSRfibers+1:end,:)))/...
+            (nHSRfibers*duration))])
+        
+        nCNneurons=size(CNoutput,1);
+        nHSRCNneuronss=nCNneurons/nANfiberTypes;
+        disp(['CN: ' num2str(sum(sum(CNoutput(end-nHSRCNneuronss+1:end,:)))...
+            /(nHSRCNneuronss*duration))])
+        disp(['IC: ' num2str(sum(sum(ICoutput))/duration)])
+        %         disp(['IC by type: ' num2str(mean(ICfiberTypeRates,2)')])
+    else
+        disp(['AN: ' num2str(mean(mean(ANprobRateOutput)))])
+        [PSTH pointsPerBin]= UTIL_makePSTH(ANprobRateOutput, dt, 0.001);
+        disp(['max max AN: ' num2str(max(max(...
+          PSTH/pointsPerBin )))])
+    end
+end
+
+
+%% figure (99) summarises main model output
+if options.showModelOutput
+    plotInstructions.figureNo=99;
+    signalRMS=mean(savedInputSignal.^2)^0.5;
+    signalRMSdb=20*log10(signalRMS/20e-6);
+    
+    % plot signal (1)
+    plotInstructions.displaydt=dt;
+    plotInstructions.numPlots=6;
+    plotInstructions.subPlotNo=1;
+    plotInstructions.title=...
+        ['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL'];
+    r=size(savedInputSignal,1);
+    if r==1, savedInputSignal=savedInputSignal'; end
+    UTIL_plotMatrix(savedInputSignal', plotInstructions);
+    
+    % stapes (2)
+    plotInstructions.subPlotNo=2;
+    plotInstructions.title= ['stapes displacement'];
+    UTIL_plotMatrix(OMEoutput, plotInstructions);
+    
+    % DRNL (3)
+    plotInstructions.subPlotNo=3;
+    plotInstructions.title= ['BM displacement'];
+    plotInstructions.yValues= savedBFlist;
+    UTIL_plotMatrix(DRNLoutput, plotInstructions);
+    
+    switch saveAN_spikesOrProbability
+        case 'spikes'
+            % AN (4)
+            plotInstructions.displaydt=ANdt;
+            plotInstructions.title='AN';
+            plotInstructions.subPlotNo=4;
+            plotInstructions.yLabel='BF';
+            plotInstructions.yValues= savedBFlist;
+            plotInstructions.rasterDotSize=1;
+            plotInstructions.plotDivider=1;
+            if sum(sum(ANoutput))<100
+                plotInstructions.rasterDotSize=3;
+            end
+            UTIL_plotMatrix(ANoutput, plotInstructions);
+            
+            % CN (5)
+            plotInstructions.displaydt=ANdt;
+            plotInstructions.subPlotNo=5;
+            plotInstructions.title='CN spikes';
+            if sum(sum(CNoutput))<100
+                plotInstructions.rasterDotSize=3;
+            end
+            UTIL_plotMatrix(CNoutput, plotInstructions);
+            
+            % IC (6)
+            plotInstructions.displaydt=ANdt;
+            plotInstructions.subPlotNo=6;
+            plotInstructions.title='IC';
+            if size(ICoutput,1)>3
+                if sum(sum(ICoutput))<100
+                    plotInstructions.rasterDotSize=3;
+                end
+                UTIL_plotMatrix(ICoutput, plotInstructions);
+            else
+                plotInstructions.title='IC (HSR) membrane potential';
+                plotInstructions.displaydt=dt;
+                plotInstructions.yLabel='V';
+                plotInstructions.zValuesRange= [-.1 0];
+                UTIL_plotMatrix(ICmembraneOutput, plotInstructions);
+            end
+            
+        otherwise % probability (4-6)
+            plotInstructions.displaydt=dt;
+            plotInstructions.numPlots=2;
+            plotInstructions.subPlotNo=2;
+            plotInstructions.yLabel='BF';
+            if nANfiberTypes>1,
+                plotInstructions.yLabel='LSR    HSR';
+                plotInstructions.plotDivider=1;
+            end
+            plotInstructions.title='AN - spike probability';
+            UTIL_plotMatrix(ANprobRateOutput, plotInstructions);
+    end
+end
+
+%% surface plot of probability
+if options.surfProbability
+    figure(97), clf
+    % select only HSR fibers at the bottom of the matrix
+    ANprobRateOutput= ANprobRateOutput(end-length(savedBFlist)+1:end,:);
+    [nY nX]=size(ANprobRateOutput);
+    if nY>2
+        time=dt*(1:nX);
+        surf(time, savedBFlist, ANprobRateOutput)
+        shading interp
+        set(gca, 'yScale','log')
+        xlim([0 max(time)]), ylim([0 max(savedBFlist)]), zlim([0 1000])
+        xlabel('time (s)')
+        ylabel('best frequency (Hz)')
+        zlabel('spike rate')
+        view([-20 60])
+        if isfield(options, 'fileName')
+            title ([options.fileName ':   ' num2str(signalRMSdb,'% 3.0f') ' dB'])
+        else
+            title ([ num2str(signalRMSdb,'% 3.0f') ' dB'])
+        end
+        
+    end
+end
+
+
+%% plot efferent control values as dB
+if options.showEfferent
+    plotInstructions=[];
+    plotInstructions.figureNo=98;
+    figure(98), clf
+    plotInstructions.displaydt=dt;
+    plotInstructions.numPlots=2;
+    plotInstructions.subPlotNo=1;
+    plotInstructions.zValuesRange=[ -25 0];
+    plotInstructions.title= ['AR strength.  Signal level= ' ...
+        num2str(signalRMSdb,'%4.0f') ' dB SPL'];
+    UTIL_plotMatrix(20*log10(ARattenuation), plotInstructions);
+    
+    plotInstructions.subPlotNo=2;
+    plotInstructions.yValues= savedBFlist;
+    plotInstructions.yLabel= 'BF';
+    plotInstructions.title= ['MOC strength'];
+    plotInstructions.zValuesRange=[ -25 0];
+    subplot(2,1,2)
+    % imagesc(MOCattenuation)
+    UTIL_plotMatrix(20*log10(MOCattenuation), plotInstructions);
+    colorbar
+end
+<<<<<<< HEAD
+
+if options.surfProbability
+%% surface plot of probability
+    figure(97), clf
+    % select only HSR fibers at the bottom of the matrix
+    ANprobRateOutput= ANprobRateOutput(end-length(savedBFlist)+1:end,:);
+    PSTHbinWidth=0.001;
+    PSTH=UTIL_PSTHmakerb(ANprobRateOutput, ANdt, PSTHbinWidth);
+    [nY nX]=size(PSTH);
+    PSTH=PSTH(:,200:end);
+    [nY nX]=size(PSTH);
+    time=PSTHbinWidth*(1:nX);
+    surf(time, savedBFlist, PSTH)
+    shading interp
+    set(gca, 'yScale','log')
+    xlim([0 max(time)])
+    ylim([0 max(savedBFlist)])
+    zlim([0 1000])
+    xlabel('time (s)')
+    ylabel('best frequency (Hz)')
+    zlabel('spike rate')
+    view([-20 60])
+%     view([0 90])
+    title ([options.fileName ':   ' num2str(signalRMSdb,'% 3.0f') ' dB'])
+end
+
+
+%% ACF plot if required
+if options.showACF
+    tic
+    method.dt=dt;
+    method.segmentNo=1;
+    method.nonlinCF=savedBFlist;
+
+    minPitch=	80; maxPitch=	4000; numPitches=100;    % specify lags
+    pitches=10.^ linspace(log10(minPitch), log10(maxPitch),numPitches);
+    pitches=fliplr(pitches);
+    filteredSACFParams.lags=1./pitches;     % autocorrelation lags vector
+    filteredSACFParams.acfTau=	.003;       % time constant of running ACF
+    filteredSACFParams.lambda=	0.12;       % slower filter to smooth ACF
+    filteredSACFParams.lambda=	0.01;       % slower filter to smooth ACF
+
+    filteredSACFParams.plotACFs=0;          % special plot (see code)
+    filteredSACFParams.plotFilteredSACF=0;  % 0 plots unfiltered ACFs
+    filteredSACFParams.plotMoviePauses=.3;          % special plot (see code)
+
+    filteredSACFParams.usePressnitzer=0; % attenuates ACF at  long lags
+    filteredSACFParams.lagsProcedure=  'useAllLags';
+    % filteredSACFParams.lagsProcedure=  'useBernsteinLagWeights';
+    % filteredSACFParams.lagsProcedure=  'omitShortLags';
+    filteredSACFParams.criterionForOmittingLags=3;
+    filteredSACFParams.plotACFsInterval=200;
+
+    if filteredSACFParams.plotACFs
+        % plot original waveform on ACF plot
+        figure(13), clf
+        subplot(4,1,1)
+=======
+    
+    %% ACF plot if required
+    if options.showACF
+        tic
+        method.dt=dt;
+        method.segmentNo=1;
+        method.nonlinCF=savedBFlist;
+        
+        minPitch=	80; maxPitch=	4000; numPitches=100;    % specify lags
+        pitches=10.^ linspace(log10(minPitch), log10(maxPitch),numPitches);
+        pitches=fliplr(pitches);
+        filteredSACFParams.lags=1./pitches;     % autocorrelation lags vector
+        filteredSACFParams.acfTau=	.003;       % time constant of running ACF
+        filteredSACFParams.lambda=	0.12;       % slower filter to smooth ACF
+        filteredSACFParams.lambda=	0.01;       % slower filter to smooth ACF
+        
+        filteredSACFParams.plotACFs=0;          % special plot (see code)
+        filteredSACFParams.plotFilteredSACF=0;  % 0 plots unfiltered ACFs
+        filteredSACFParams.plotMoviePauses=.3;          % special plot (see code)
+        
+        filteredSACFParams.usePressnitzer=0; % attenuates ACF at  long lags
+        filteredSACFParams.lagsProcedure=  'useAllLags';
+        % filteredSACFParams.lagsProcedure=  'useBernsteinLagWeights';
+        % filteredSACFParams.lagsProcedure=  'omitShortLags';
+        filteredSACFParams.criterionForOmittingLags=3;
+        filteredSACFParams.plotACFsInterval=200;
+        
+        if filteredSACFParams.plotACFs
+            % plot original waveform on ACF plot
+            figure(13), clf
+            subplot(4,1,1)
+            t=dt*(1:length(savedInputSignal));
+            plot(t,savedInputSignal)
+            xlim([0 t(end)])
+            title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']);
+        end
+        
+        % plot original waveform on summary/smoothed ACF plot
+        figure(96), clf
+        subplot(2,1,1)
+>>>>>>> 77bdf847c4be95767d3ea519d77878f9414a2686
+        t=dt*(1:length(savedInputSignal));
+        plot(t,savedInputSignal)
+        xlim([0 t(end)])
+        title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']);
+        
+        
+        % compute ACF
+        switch saveAN_spikesOrProbability
+            case 'probability'
+                inputToACF=ANprobRateOutput.^0.5;
+            otherwise
+                inputToACF=ANoutput;
+        end
+        
+        disp ('computing ACF...')
+        [P, BFlist, sacf, boundaryValue] = ...
+            filteredSACF(inputToACF, method, filteredSACFParams);
+        disp(' ACF done.')
+        
+        % SACF
+        subplot(2,1,2)
+        imagesc(P)
+        ylabel('periodicities (Hz)')
+        xlabel('time (s)')
+        title(['running smoothed (root) SACF. ' saveAN_spikesOrProbability ' input'])
+        pt=[1 get(gca,'ytick')]; % force top xtick to show
+        set(gca,'ytick',pt)
+        set(gca,'ytickLabel', round(pitches(pt)))
+        tt=get(gca,'xtick');
+        set(gca,'xtickLabel', round(100*t(tt))/100)
+    end
+    
+    path(restorePath)
diff -r fafe69c43108 -r 45f28c49461e old files/test_MAP1_14.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/old files/test_MAP1_14.m	Mon Jun 13 18:21:05 2011 +0100
@@ -0,0 +1,222 @@
+function test_MAP1_14
+% test_MAP1_14 is a general purpose test routine that can be adjusted to
+% test a number of different applications of MAP1_14
+%
+% A range of options are supplied in the early part of the program
+%
+% One use of the function is to create demonstrations; filenames <demoxx>
+%  to illustrate particular features
+%
+% #1
+% Identify the file (in 'MAPparamsName') containing the model parameters
+%
+% #2
+% Identify the kind of model required (in 'AN_spikesOrProbability').
+%  A full brainstem model (spikes) can be computed or a shorter model
+%  (probability) that computes only so far as the auditory nerve
+%
+% #3
+% Choose between a tone signal or file input (in 'signalType')
+%
+% #4
+% Set the signal rms level (in leveldBSPL)
+%
+% #5
+% Indentify the channels in terms of their best frequencies in the vector
+%  BFlist.
+%
+% Last minute changes to the parameters fetched earlier can be made using
+%  the cell array of strings 'paramChanges'.
+%  Each string must have the same format as the corresponding line in the
+%  file identified in 'MAPparamsName'
+%
+% When the demonstration is satisfactory, freeze it by renaming it <demoxx>
+
+
+%%  #1 parameter file name
+MAPparamsName='Normal';
+
+
+%% #2 probability (fast) or spikes (slow) representation
+AN_spikesOrProbability='spikes';
+% or
+AN_spikesOrProbability='probability';
+
+
+%% #3 pure tone, harmonic sequence or speech file input
+signalType= 'tones';
+duration=0.100;                 % seconds
+% duration=0.020;                 % seconds
+sampleRate= 100000;
+% toneFrequency= 250:250:8000;    % harmonic sequence (Hz)
+toneFrequency= 1000;            % or a pure tone (Hz8
+
+rampDuration=.005;              % seconds
+
+% % or
+% signalType= 'file';
+% fileName='twister_44kHz';
+% % fileName='new-da-44khz';
+
+% ? and mix with an optional second file?
+mixerFile=[];
+%or
+<<<<<<< HEAD
+% mixerFile='babble';
+% leveldBSPL2=30;
+=======
+mixerFile='babble';
+leveldBSPL2=60;
+>>>>>>> 77bdf847c4be95767d3ea519d77878f9414a2686
+
+%% #4 rms level
+% signal details
+leveldBSPL= 70;                  % dB SPL
+
+
+%% #5 number of channels in the model
+%   21-channel model (log spacing)
+numChannels=21;
+lowestBF=250; 	highestBF= 8000;
+BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels));
+
+%   or specify your own channel BFs
+% BFlist=toneFrequency;
+
+
+%% #6 change model parameters
+paramChanges=[];
+
+% or
+% Parameter changes can be used to change one or more model parameters
+%  *after* the MAPparams file has been read
+% This example declares only one fiber type with a calcium clearance time
+% constant of 80e-6 s (HSR fiber) when the probability option is selected.
+% paramChanges={'AN_IHCsynapseParams.ANspeedUpFactor=5;', ...
+%     'IHCpreSynapseParams.tauCa=86e-6;'};
+% paramChanges={'DRNLParams.rateToAttenuationFactorProb = 0;'};
+
+%% delare showMap options
+global showMapOptions
+
+% or (example: show everything including an smoothed SACF output
+showMapOptions.showModelParameters=1;
+showMapOptions.showModelOutput=1;
+showMapOptions.printFiringRates=1;
+showMapOptions.showACF=0;
+showMapOptions.showEfferent=1;
+if strcmp(AN_spikesOrProbability, 'probability')
+    showMapOptions.surfProbability=1;
+else
+    showMapOptions.surfProbability=0;
+end
+if strcmp(signalType, 'file')
+    showMapOptions.fileName=fileName;
+else
+    showMapOptions.fileName=[];
+end
+
+%% Generate stimuli
+
+dbstop if error
+restorePath=path;
+addpath (['..' filesep 'MAP'],    ['..' filesep 'wavFileStore'])
+switch signalType
+    case 'tones'
+        inputSignal=createMultiTone(sampleRate, toneFrequency, ...
+            leveldBSPL, duration, rampDuration);
+        
+    case 'file'
+        %% file input simple or mixed
+        [inputSignal sampleRate]=wavread(fileName);
+            dt=1/sampleRate;
+        inputSignal=inputSignal(:,1);
+        targetRMS=20e-6*10^(leveldBSPL/20);
+        rms=(mean(inputSignal.^2))^0.5;
+        amp=targetRMS/rms;
+        inputSignal=inputSignal*amp;
+<<<<<<< HEAD
+        silence=zeros(10000,1);
+        inputSignal=[silence; inputSignal];
+=======
+        silence= zeros(1,round(0.1/dt));
+        inputSignal= [silence inputSignal' silence];
+        
+>>>>>>> 77bdf847c4be95767d3ea519d77878f9414a2686
+        if ~isempty(mixerFile)
+            [inputSignal2 sampleRate2]=wavread(mixerFile);
+            if ~isequal(sampleRate,sampleRate2)
+                error...
+                    ('file and mixer file have different sample rates')
+            end
+            inputSignal2=inputSignal2(:,1);
+            [r c]=size(inputSignal);
+            inputSignal2=inputSignal2(1:r);
+            targetRMS=20e-6*10^(leveldBSPL2/20);
+            rms=(mean(inputSignal2.^2))^0.5;
+            amp=targetRMS/rms;
+            inputSignal2=inputSignal2*amp;
+            rampDuration=dt*length(inputSignal2)/3;
+            if rampDuration>0.5*duration, rampDuration=duration/2; end
+            rampTime=dt:dt:rampDuration;
+            time=dt*(1:length(inputSignal2));
+            ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
+                ones(1,length(time)-length(rampTime))];
+            inputSignal2=inputSignal2.*ramp';
+            ramp=fliplr(ramp);
+            inputSignal2=inputSignal2.*ramp';
+        end
+        
+        inputSignal=inputSignal+inputSignal2;
+end
+
+
+%% run the model
+tic
+
+fprintf('\n')
+disp(['Signal duration= ' num2str(length(inputSignal)/sampleRate)])
+disp([num2str(numChannels) ' channel model'])
+disp('Computing ...')
+
+restorePath=path;
+addpath (['..' filesep 'MAP'])
+
+MAP1_14(inputSignal, sampleRate, BFlist, ...
+    MAPparamsName, AN_spikesOrProbability, paramChanges);
+path(restorePath)
+toc
+
+% the model run is now complete. Now display the results
+showMAP
+for i=1:length(paramChanges)
+    disp(paramChanges{i})
+end
+
+toc
+path(restorePath)
+
+
+function inputSignal=createMultiTone(sampleRate, toneFrequency, ...
+    leveldBSPL, duration, rampDuration)
+% Create pure tone stimulus
+dt=1/sampleRate; % seconds
+time=dt: dt: duration;
+inputSignal=sum(sin(2*pi*toneFrequency'*time), 1);
+amp=10^(leveldBSPL/20)*28e-6;   % converts to Pascals (peak)
+inputSignal=amp*inputSignal;
+
+% apply ramps
+% catch rampTime error
+if rampDuration>0.5*duration, rampDuration=duration/2; end
+rampTime=dt:dt:rampDuration;
+ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
+    ones(1,length(time)-length(rampTime))];
+inputSignal=inputSignal.*ramp;
+ramp=fliplr(ramp);
+inputSignal=inputSignal.*ramp;
+
+% add 10 ms silence
+silence= zeros(1,round(0.03/dt));
+inputSignal= [silence inputSignal silence];
+
diff -r fafe69c43108 -r 45f28c49461e parameterStore/MAPparamsNormal.m
--- a/parameterStore/MAPparamsNormal.m	Mon Jun 13 17:30:57 2011 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,339 +0,0 @@
-function method=MAPparamsNormal ...
-    (BFlist, sampleRate, showParams)
-% MAPparams<> establishes a complete set of MAP parameters
-% Parameter file names must be of the form <MAPparams> <name>
-%
-% input arguments
-%  BFlist     (optional) specifies the desired list of channel BFs
-%    otherwise defaults set below
-%  sampleRate (optional), default is 50000.
-%  showParams (optional) =1 prints out the complete set of parameters
-% output argument
-%  method passes a miscelleny of values
-
-global inputStimulusParams OMEParams DRNLParams
-global IHC_VResp_VivoParams IHCpreSynapseParams  AN_IHCsynapseParams
-global MacGregorParams MacGregorMultiParams  filteredSACFParams
-global experiment % used by calls from multiThreshold only
-global IHC_cilia_RPParams
-
-currentFile=mfilename;                      % i.e. the name of this mfile
-method.parameterSource=currentFile(10:end); % for the record
-
-efferentDelay=0.010;
-method.segmentDuration=efferentDelay;
-
-if nargin<3, showParams=0; end
-if nargin<2, sampleRate=50000; end
-if nargin<1 || BFlist(1)<0 % if BFlist= -1, set BFlist to default
-    lowestBF=250; 	highestBF= 8000; 	numChannels=21;
-    % 21 chs (250-8k)includes BFs at 250 500 1000 2000 4000 8000
-    BFlist=round(logspace(log10(lowestBF),log10(highestBF),numChannels));
-end
-% BFlist=1000;
-
-% preserve for backward campatibility
-method.nonlinCF=BFlist; 
-method.dt=1/sampleRate; 
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% set  model parameters
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-%%  #1 inputStimulus
-inputStimulusParams=[];
-inputStimulusParams.sampleRate= sampleRate; 
-
-%%  #2 outerMiddleEar
-OMEParams=[];  % clear the structure first
-% outer ear resonances band pass filter  [gain lp order  hp]
-OMEParams.externalResonanceFilters=      [ 10 1 1000 4000];
-                                       
-% highpass stapes filter  
-%  Huber gives 2e-9 m at 80 dB and 1 kHz (2e-13 at 0 dB SPL)
-OMEParams.OMEstapesLPcutoff= 1000;
-OMEParams.stapesScalar=	     45e-9;
-
-% 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.004;   % * N(all ICspikes)
-%     OMEParams.rateToAttenuationFactor=0;   % * N(all ICspikes)
-
-% 'probability model': Ar based on AN firing probabilities (LSR)
-OMEParams.rateToAttenuationFactorProb=0.003;% * 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;
-
-%%  #3 DRNL
-DRNLParams=[];  % clear the structure first
-DRNLParams.BFlist=BFlist;
-
-% DRNL nonlinear path
-DRNLParams.a=3e4;     % nonlinear path gain (below compression threshold)
-DRNLParams.a=5e2;     % 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.c=0.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)
-DRNLParams.nlBWs=  p * BFlist + q;
-DRNLParams.p=p;   DRNLParams.q=q;   % save p and q for printing only
-
-% DRNL linear path:
-DRNLParams.g=100;     % linear path gain factor
-% linCF is not necessarily the same as nonlinCF
-minLinCF=153.13; coeffLinCF=0.7341;   % linCF>nonlinBF for BF < 1 kHz
-DRNLParams.linCFs=minLinCF+coeffLinCF*BFlist;
-DRNLParams.linOrder=	3;           % order of linear gammatone filters
-minLinBW=100; coeffLinBW=0.6531;
-DRNLParams.linBWs=minLinBW + coeffLinBW*BFlist; % bandwidths of linear  filters
-
-% DRNL MOC efferents
-DRNLParams.MOCdelay = efferentDelay;            % must be < segment length!
-% 'spikes' model: MOC based on brainstem spiking activity (HSR)
-DRNLParams.rateToAttenuationFactor = .009;  % strength of MOC
-DRNLParams.rateToAttenuationFactor = .004;  % strength of MOC
-%      DRNLParams.rateToAttenuationFactor = 0;  % strength of MOC
-
-% 'probability' model: MOC based on AN spiking activity (HSR)
-DRNLParams.rateToAttenuationFactorProb = .004;  % strength of MOC
-% DRNLParams.rateToAttenuationFactorProb = .0;  % strength of MOC
-DRNLParams.MOCtau =.1;                         % smoothing for MOC
-DRNLParams.MOCrateThreshold =50;                % set to AN rate threshold
-
-
-%% #4 IHC_cilia_RPParams
-
-IHC_cilia_RPParams.tc=	0.0003;   % 0.0003 filter time simulates viscocity
-% IHC_cilia_RPParams.tc=	0.0005;   % 0.0003 filter time simulates viscocity
-IHC_cilia_RPParams.C=	0.05;      % 0.1 scalar (C_cilia ) 
-IHC_cilia_RPParams.u0=	5e-9;       
-IHC_cilia_RPParams.s0=	30e-9;
-IHC_cilia_RPParams.u1=	1e-9;
-IHC_cilia_RPParams.s1=	1e-9;
-
-IHC_cilia_RPParams.Gmax= 5e-9;    % 2.5e-9 maximum conductance (Siemens)
-IHC_cilia_RPParams.Ga=	1e-9;  % 4.3e-9 fixed apical membrane conductance
-
-%  #5 IHC_RP
-IHC_cilia_RPParams.Cab=	4e-012;         % IHC capacitance (F)
-IHC_cilia_RPParams.Cab=	1e-012;         % IHC capacitance (F)
-IHC_cilia_RPParams.Et=	0.100;          % endocochlear potential (V)
-
-IHC_cilia_RPParams.Gk=	2e-008;         % 1e-8 potassium conductance (S)
-IHC_cilia_RPParams.Ek=	-0.08;          % -0.084 K equilibrium potential
-IHC_cilia_RPParams.Rpc=	0.04;           % combined resistances
-
-
-%%  #5 IHCpreSynapse
-IHCpreSynapseParams=[];
-IHCpreSynapseParams.GmaxCa=	14e-9;% maximum calcium conductance
-IHCpreSynapseParams.GmaxCa=	12e-9;% maximum calcium conductance
-IHCpreSynapseParams.ECa=	0.066;  % calcium equilibrium potential
-IHCpreSynapseParams.beta=	400;	% determine Ca channel opening
-IHCpreSynapseParams.gamma=	100;	% determine Ca channel opening
-IHCpreSynapseParams.tauM=	0.00005;	% membrane time constant ?0.1ms
-IHCpreSynapseParams.power=	3;
-% 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=70e-6;            % seconds
-IHCpreSynapseParams.tauCa= [LSRtauCa HSRtauCa]; %LSR and HSR fiber
-
-%%  #6 AN_IHCsynapse
-% c=kym/(y(l+r)+kl)	(spontaneous rate)
-% c=(approx)  ym/l  (saturated rate)
-AN_IHCsynapseParams=[];             % clear the structure first
-AN_IHCsynapseParams.M=	12;         % maximum vesicles at synapse
-AN_IHCsynapseParams.y=	4;          % depleted vesicle replacement rate
-AN_IHCsynapseParams.y=	6;          % depleted vesicle replacement rate
-
-AN_IHCsynapseParams.x=	30;         % replenishment from re-uptake store
-AN_IHCsynapseParams.x=	60;         % replenishment from re-uptake store
-
-% reduce l to increase saturated rate
-AN_IHCsynapseParams.l=	100; % *loss rate of vesicles from the cleft
-AN_IHCsynapseParams.l=	250; % *loss rate of vesicles from the cleft
-
-AN_IHCsynapseParams.r=	500; % *reuptake rate from cleft into cell
-% AN_IHCsynapseParams.r=	300; % *reuptake rate from cleft into cell
-
-AN_IHCsynapseParams.refractory_period=	0.00075;
-% number of AN fibers at each BF (used only for spike generation)
-AN_IHCsynapseParams.numFibers=	100; 
-AN_IHCsynapseParams.TWdelay=0.004;  % ?delay before stimulus first spike
-
-AN_IHCsynapseParams.ANspeedUpFactor=5; % longer epochs for computing spikes.
-
-%%  #7 MacGregorMulti (first order brainstem neurons)
-MacGregorMultiParams=[];
-MacGregorMultiType='chopper'; % MacGregorMultiType='primary-like'; %choose
-switch MacGregorMultiType
-    case 'primary-like'
-        MacGregorMultiParams.nNeuronsPerBF=	10;   % N neurons per BF
-        MacGregorMultiParams.type = 'primary-like cell';
-        MacGregorMultiParams.fibersPerNeuron=4;   % N input fibers
-        MacGregorMultiParams.dendriteLPfreq=200;  % dendritic filter
-        MacGregorMultiParams.currentPerSpike=0.11e-6; % (A) per spike
-        MacGregorMultiParams.Cap=4.55e-9;   % cell capacitance (Siemens)
-        MacGregorMultiParams.tauM=5e-4;     % membrane time constant (s)
-        MacGregorMultiParams.Ek=-0.01;      % K+ eq. potential (V)
-        MacGregorMultiParams.dGkSpike=3.64e-5; % K+ cond.shift on spike,S
-        MacGregorMultiParams.tauGk=	0.0012; % K+ conductance tau (s)
-        MacGregorMultiParams.Th0=	0.01;   % equilibrium threshold (V)
-        MacGregorMultiParams.c=	0.01;       % threshold shift on spike, (V)
-        MacGregorMultiParams.tauTh=	0.015;  % variable threshold tau
-        MacGregorMultiParams.Er=-0.06;      % resting potential (V)
-        MacGregorMultiParams.Eb=0.06;       % spike height (V)
-
-    case 'chopper'
-        MacGregorMultiParams.nNeuronsPerBF=	10;   % N neurons per BF
-        MacGregorMultiParams.type = 'chopper cell';
-        MacGregorMultiParams.fibersPerNeuron=10;  % N input fibers
-%         MacGregorMultiParams.fibersPerNeuron=6;  % N input fibers
-
-        MacGregorMultiParams.dendriteLPfreq=50;   % dendritic filter
-        MacGregorMultiParams.currentPerSpike=35e-9; % *per spike
-        MacGregorMultiParams.currentPerSpike=30e-9; % *per spike
-        
-        MacGregorMultiParams.Cap=1.67e-8; % ??cell capacitance (Siemens)
-        MacGregorMultiParams.tauM=0.002;  % membrane time constant (s)
-        MacGregorMultiParams.Ek=-0.01;    % K+ eq. potential (V)
-        MacGregorMultiParams.dGkSpike=1.33e-4; % K+ cond.shift on spike,S
-        MacGregorMultiParams.tauGk=	0.0005;% K+ conductance tau (s)
-        MacGregorMultiParams.Th0=	0.01; % equilibrium threshold (V)
-        MacGregorMultiParams.c=	0;        % threshold shift on spike, (V)
-        MacGregorMultiParams.tauTh=	0.02; % variable threshold tau
-        MacGregorMultiParams.Er=-0.06;    % resting potential (V)
-        MacGregorMultiParams.Eb=0.06;     % spike height (V)
-        MacGregorMultiParams.PSTHbinWidth=	1e-4;
-end
-
-%%  #8 MacGregor (second-order neuron). Only one per channel
-MacGregorParams=[];                 % clear the structure first
-MacGregorParams.type = 'chopper cell';
-MacGregorParams.fibersPerNeuron=10; % N input fibers
-MacGregorParams.dendriteLPfreq=100; % dendritic filter
-MacGregorParams.currentPerSpike=120e-9;% *(A) per spike
-MacGregorParams.currentPerSpike=30e-9;% *(A) per spike
-
-MacGregorParams.Cap=16.7e-9;        % cell capacitance (Siemens)
-MacGregorParams.tauM=0.002;         % membrane time constant (s)
-MacGregorParams.Ek=-0.01;           % K+ eq. potential (V)
-MacGregorParams.dGkSpike=1.33e-4;   % K+ cond.shift on spike,S
-MacGregorParams.tauGk=	0.0005;     % K+ conductance tau (s)
-MacGregorParams.Th0=	0.01;       % equilibrium threshold (V)
-MacGregorParams.c=	0;              % threshold shift on spike, (V)
-MacGregorParams.tauTh=	0.02;       % variable threshold tau
-MacGregorParams.Er=-0.06;           % resting potential (V)
-MacGregorParams.Eb=0.06;            % spike height (V)
-MacGregorParams.debugging=0;        % (special)
-% wideband accepts input from all channels (of same fiber type)
-% use wideband to create inhibitory units
-MacGregorParams.wideband=0;         % special for wideband units
-% MacGregorParams.saveAllData=0;
-
-%%  #9 filteredSACF
-minPitch=	300; maxPitch=	3000; numPitches=60;    % specify lags
-pitches=100*log10(logspace(minPitch/100, maxPitch/100, numPitches));
-filteredSACFParams.lags=1./pitches;     % autocorrelation lags vector
-filteredSACFParams.acfTau=	.003;       % time constant of running ACF
-filteredSACFParams.lambda=	0.12;       % slower filter to smooth ACF
-filteredSACFParams.plotFilteredSACF=1;  % 0 plots unfiltered ACFs
-filteredSACFParams.plotACFs=0;          % special plot (see code)
-%  filteredSACFParams.usePressnitzer=0; % attenuates ACF at  long lags
-filteredSACFParams.lagsProcedure=  'useAllLags';
-% filteredSACFParams.lagsProcedure=  'useBernsteinLagWeights';
-% filteredSACFParams.lagsProcedure=  'omitShortLags';
-filteredSACFParams.criterionForOmittingLags=3;
-
-% checks
-if AN_IHCsynapseParams.numFibers<MacGregorMultiParams.fibersPerNeuron
-    error('MacGregorMulti: too few input fibers for input to MacG unit')
-end
-
-
-%% write all parameters to the command window
-% showParams is currently set at the top of htis function
-if showParams
-    fprintf('\n %%%%%%%%\n')
-    fprintf('\n%s\n', method.parameterSource)
-    fprintf('\n')
-    nm=UTIL_paramsList(whos);
-    for i=1:length(nm)
-        %         eval(['UTIL_showStruct(' nm{i} ', ''' nm{i} ''')'])
-        if ~strcmp(nm(i), 'method')
-            eval(['UTIL_showStructureSummary(' nm{i} ', ''' nm{i} ''', 10)'])
-        end
-    end
-end
-
-
-
-% **********************************************************************  comparison data
-% store individual data here for display on the multiThreshold GUI (if used)
-% the final value in each vector is an identifier (BF or duration))
-if isstruct(experiment)
-    switch experiment.paradigm
-        case {'IFMC','IFMC_8ms'}
-            % based on MPa
-            comparisonData=[
-                66	51	49	48	46	45	54	250;
-                60	54	46	42	39	49	65	500;
-                64	51	38	32	33	59	75	1000;
-                59	51	36	30	41	81	93	2000;
-                71	63	53	44	36	76	95	4000;
-                70	64	43	35	35	66	88	6000;
-                110	110	110	110	110	110	110	8000;
-                ];
-            if length(BFlist)==1 && ~isempty(comparisonData)
-                availableFrequencies=comparisonData(:,end)';
-                findRow= find(BFlist==availableFrequencies);
-                if ~isempty (findRow)
-                    experiment.comparisonData=comparisonData(findRow,:);
-                end
-            end
-
-        case {'TMC','TMC_8ms'}
-            % based on MPa
-            comparisonData=[
-                48	58	63	68	75	80	85	92	99	250;
-                33	39	40	49	52	61	64	77	79	500;
-                39	42	50	81	83	92	96	97	110	1000;
-                24	26	32	37	46	51	59	71	78	2000;
-                65	68	77	85	91	93	110	110	110	4000;
-                20	19	26	44	80	95	96	110	110	6000;
-                ];
-            if length(BFlist)==1 && ~isempty(comparisonData)
-                availableFrequencies=comparisonData(:,end)';
-                findRow= find(BFlist==availableFrequencies);
-                if ~isempty (findRow)
-                    experiment.comparisonData=comparisonData(findRow,:);
-                end
-            end
-
-        case { 'absThreshold',  'absThreshold_8'}
-            % MPa thresholds
-            experiment.comparisonData=[
-                32	26	16	18	22	22 0.008;
-                16	13	6	9	15	11 0.500
-                ];
-
-
-        otherwise
-            experiment.comparisonData=[];
-    end
-end
-
-
diff -r fafe69c43108 -r 45f28c49461e testPrograms/showMAP.m
--- a/testPrograms/showMAP.m	Mon Jun 13 17:30:57 2011 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,275 +0,0 @@
-function showMAP (options)
-% defaults
-% options.showModelParameters=1;
-% options.showModelOutput=1;
-% options.printFiringRates=1;
-% options.showACF=1;
-% options.showEfferent=1;
-% options.surfProbability=0;
-% options.fileName=[];
-
-dbstop if warning
-
-global dt ANdt saveAN_spikesOrProbability savedBFlist saveMAPparamsName...
-    savedInputSignal TMoutput OMEoutput ARattenuation ...
-    DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
-    IHCoutput ANprobRateOutput ANoutput savePavailable tauCas  ...
-    CNoutput  ICoutput ICmembraneOutput ICfiberTypeRates MOCattenuation
-global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams
-global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams
-
-
-restorePath=path;
-addpath ( ['..' filesep 'utilities'], ['..' filesep 'parameterStore'])
-
-if nargin<1
-    options.showModelParameters=1;
-    options.showModelOutput=1;
-    options.printFiringRates=1;
-    options.showACF=0;
-    options.showEfferent=1;
-    options.surfProbability=0;
-    options.fileName=[];
-end
-
-if options.showModelParameters
-    % Read parameters from MAPparams<***> file in 'parameterStore' folder
-    %  and print out all parameters
-    cmd=['MAPparams' saveMAPparamsName ...
-        '(-1, 1/dt, 1);'];
-    eval(cmd);
-end
-
-if options.printFiringRates
-    %% print summary firing rates
-    fprintf('\n\n')
-    disp('summary')
-    disp(['AR: ' num2str(min(ARattenuation))])
-    disp(['MOC: ' num2str(min(min(MOCattenuation)))])
-    nANfiberTypes=length(tauCas);
-    if strcmp(saveAN_spikesOrProbability, 'spikes')
-        nANfibers=size(ANoutput,1);
-        nHSRfibers=nANfibers/nANfiberTypes;
-        duration=size(TMoutput,2)*dt;
-        disp(['AN: ' num2str(sum(sum(ANoutput(end-nHSRfibers+1:end,:)))/...
-            (nHSRfibers*duration))])
-        
-        nCNneurons=size(CNoutput,1);
-        nHSRCNneuronss=nCNneurons/nANfiberTypes;
-        disp(['CN: ' num2str(sum(sum(CNoutput(end-nHSRCNneuronss+1:end,:)))...
-            /(nHSRCNneuronss*duration))])
-        disp(['IC: ' num2str(sum(sum(ICoutput))/duration)])
-        %         disp(['IC by type: ' num2str(mean(ICfiberTypeRates,2)')])
-    else
-        disp(['AN: ' num2str(mean(mean(ANprobRateOutput)))])
-        [PSTH pointsPerBin]= UTIL_makePSTH(ANprobRateOutput, dt, 0.001);
-        disp(['max max AN: ' num2str(max(max(...
-          PSTH/pointsPerBin )))])
-    end
-end
-
-
-%% figure (99) summarises main model output
-if options.showModelOutput
-    plotInstructions.figureNo=99;
-    signalRMS=mean(savedInputSignal.^2)^0.5;
-    signalRMSdb=20*log10(signalRMS/20e-6);
-    
-    % plot signal (1)
-    plotInstructions.displaydt=dt;
-    plotInstructions.numPlots=6;
-    plotInstructions.subPlotNo=1;
-    plotInstructions.title=...
-        ['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL'];
-    r=size(savedInputSignal,1);
-    if r==1, savedInputSignal=savedInputSignal'; end
-    UTIL_plotMatrix(savedInputSignal', plotInstructions);
-    
-    % stapes (2)
-    plotInstructions.subPlotNo=2;
-    plotInstructions.title= ['stapes displacement'];
-    UTIL_plotMatrix(OMEoutput, plotInstructions);
-    
-    % DRNL (3)
-    plotInstructions.subPlotNo=3;
-    plotInstructions.title= ['BM displacement'];
-    plotInstructions.yValues= savedBFlist;
-    UTIL_plotMatrix(DRNLoutput, plotInstructions);
-    
-    switch saveAN_spikesOrProbability
-        case 'spikes'
-            % AN (4)
-            plotInstructions.displaydt=ANdt;
-            plotInstructions.title='AN';
-            plotInstructions.subPlotNo=4;
-            plotInstructions.yLabel='BF';
-            plotInstructions.yValues= savedBFlist;
-            plotInstructions.rasterDotSize=1;
-            plotInstructions.plotDivider=1;
-            if sum(sum(ANoutput))<100
-                plotInstructions.rasterDotSize=3;
-            end
-            UTIL_plotMatrix(ANoutput, plotInstructions);
-            
-            % CN (5)
-            plotInstructions.displaydt=ANdt;
-            plotInstructions.subPlotNo=5;
-            plotInstructions.title='CN spikes';
-            if sum(sum(CNoutput))<100
-                plotInstructions.rasterDotSize=3;
-            end
-            UTIL_plotMatrix(CNoutput, plotInstructions);
-            
-            % IC (6)
-            plotInstructions.displaydt=ANdt;
-            plotInstructions.subPlotNo=6;
-            plotInstructions.title='IC';
-            if size(ICoutput,1)>3
-                if sum(sum(ICoutput))<100
-                    plotInstructions.rasterDotSize=3;
-                end
-                UTIL_plotMatrix(ICoutput, plotInstructions);
-            else
-                plotInstructions.title='IC (HSR) membrane potential';
-                plotInstructions.displaydt=dt;
-                plotInstructions.yLabel='V';
-                plotInstructions.zValuesRange= [-.1 0];
-                UTIL_plotMatrix(ICmembraneOutput, plotInstructions);
-            end
-            
-        otherwise % probability (4-6)
-            plotInstructions.displaydt=dt;
-            plotInstructions.numPlots=2;
-            plotInstructions.subPlotNo=2;
-            plotInstructions.yLabel='BF';
-            if nANfiberTypes>1,
-                plotInstructions.yLabel='LSR    HSR';
-                plotInstructions.plotDivider=1;
-            end
-            plotInstructions.title='AN - spike probability';
-            UTIL_plotMatrix(ANprobRateOutput, plotInstructions);
-    end
-end
-
-%% surface plot of probability
-if options.surfProbability
-    figure(97), clf
-    % select only HSR fibers at the bottom of the matrix
-    ANprobRateOutput= ANprobRateOutput(end-length(savedBFlist)+1:end,:);
-    [nY nX]=size(ANprobRateOutput);
-    if nY>2
-        time=dt*(1:nX);
-        surf(time, savedBFlist, ANprobRateOutput)
-        shading interp
-        set(gca, 'yScale','log')
-        xlim([0 max(time)]), ylim([0 max(savedBFlist)]), zlim([0 1000])
-        xlabel('time (s)')
-        ylabel('best frequency (Hz)')
-        zlabel('spike rate')
-        view([-20 60])
-        if isfield(options, 'fileName')
-            title ([options.fileName ':   ' num2str(signalRMSdb,'% 3.0f') ' dB'])
-        else
-            title ([ num2str(signalRMSdb,'% 3.0f') ' dB'])
-        end
-        
-    end
-end
-
-
-%% plot efferent control values as dB
-if options.showEfferent
-    plotInstructions=[];
-    plotInstructions.figureNo=98;
-    figure(98), clf
-    plotInstructions.displaydt=dt;
-    plotInstructions.numPlots=2;
-    plotInstructions.subPlotNo=1;
-    plotInstructions.zValuesRange=[ -25 0];
-    plotInstructions.title= ['AR strength.  Signal level= ' ...
-        num2str(signalRMSdb,'%4.0f') ' dB SPL'];
-    UTIL_plotMatrix(20*log10(ARattenuation), plotInstructions);
-    
-    plotInstructions.subPlotNo=2;
-    plotInstructions.yValues= savedBFlist;
-    plotInstructions.yLabel= 'BF';
-    plotInstructions.title= ['MOC strength'];
-    plotInstructions.zValuesRange=[ -25 0];
-    subplot(2,1,2)
-    % imagesc(MOCattenuation)
-    UTIL_plotMatrix(20*log10(MOCattenuation), plotInstructions);
-    colorbar
-end
-    
-    %% ACF plot if required
-    if options.showACF
-        tic
-        method.dt=dt;
-        method.segmentNo=1;
-        method.nonlinCF=savedBFlist;
-        
-        minPitch=	80; maxPitch=	4000; numPitches=100;    % specify lags
-        pitches=10.^ linspace(log10(minPitch), log10(maxPitch),numPitches);
-        pitches=fliplr(pitches);
-        filteredSACFParams.lags=1./pitches;     % autocorrelation lags vector
-        filteredSACFParams.acfTau=	.003;       % time constant of running ACF
-        filteredSACFParams.lambda=	0.12;       % slower filter to smooth ACF
-        filteredSACFParams.lambda=	0.01;       % slower filter to smooth ACF
-        
-        filteredSACFParams.plotACFs=0;          % special plot (see code)
-        filteredSACFParams.plotFilteredSACF=0;  % 0 plots unfiltered ACFs
-        filteredSACFParams.plotMoviePauses=.3;          % special plot (see code)
-        
-        filteredSACFParams.usePressnitzer=0; % attenuates ACF at  long lags
-        filteredSACFParams.lagsProcedure=  'useAllLags';
-        % filteredSACFParams.lagsProcedure=  'useBernsteinLagWeights';
-        % filteredSACFParams.lagsProcedure=  'omitShortLags';
-        filteredSACFParams.criterionForOmittingLags=3;
-        filteredSACFParams.plotACFsInterval=200;
-        
-        if filteredSACFParams.plotACFs
-            % plot original waveform on ACF plot
-            figure(13), clf
-            subplot(4,1,1)
-            t=dt*(1:length(savedInputSignal));
-            plot(t,savedInputSignal)
-            xlim([0 t(end)])
-            title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']);
-        end
-        
-        % plot original waveform on summary/smoothed ACF plot
-        figure(96), clf
-        subplot(2,1,1)
-        t=dt*(1:length(savedInputSignal));
-        plot(t,savedInputSignal)
-        xlim([0 t(end)])
-        title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']);
-        
-        
-        % compute ACF
-        switch saveAN_spikesOrProbability
-            case 'probability'
-                inputToACF=ANprobRateOutput.^0.5;
-            otherwise
-                inputToACF=ANoutput;
-        end
-        
-        disp ('computing ACF...')
-        [P, BFlist, sacf, boundaryValue] = ...
-            filteredSACF(inputToACF, method, filteredSACFParams);
-        disp(' ACF done.')
-        
-        % SACF
-        subplot(2,1,2)
-        imagesc(P)
-        ylabel('periodicities (Hz)')
-        xlabel('time (s)')
-        title(['running smoothed (root) SACF. ' saveAN_spikesOrProbability ' input'])
-        pt=[1 get(gca,'ytick')]; % force top xtick to show
-        set(gca,'ytick',pt)
-        set(gca,'ytickLabel', round(pitches(pt)))
-        tt=get(gca,'xtick');
-        set(gca,'xtickLabel', round(100*t(tt))/100)
-    end
-    
-    path(restorePath)
diff -r fafe69c43108 -r 45f28c49461e testPrograms/test_MAP1_14.m
--- a/testPrograms/test_MAP1_14.m	Mon Jun 13 17:30:57 2011 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,208 +0,0 @@
-function test_MAP1_14
-% test_MAP1_14 is a general purpose test routine that can be adjusted to
-% test a number of different applications of MAP1_14
-%
-% A range of options are supplied in the early part of the program
-%
-% One use of the function is to create demonstrations; filenames <demoxx>
-%  to illustrate particular features
-%
-% #1
-% Identify the file (in 'MAPparamsName') containing the model parameters
-%
-% #2
-% Identify the kind of model required (in 'AN_spikesOrProbability').
-%  A full brainstem model (spikes) can be computed or a shorter model
-%  (probability) that computes only so far as the auditory nerve
-%
-% #3
-% Choose between a tone signal or file input (in 'signalType')
-%
-% #4
-% Set the signal rms level (in leveldBSPL)
-%
-% #5
-% Indentify the channels in terms of their best frequencies in the vector
-%  BFlist.
-%
-% Last minute changes to the parameters fetched earlier can be made using
-%  the cell array of strings 'paramChanges'.
-%  Each string must have the same format as the corresponding line in the
-%  file identified in 'MAPparamsName'
-%
-% When the demonstration is satisfactory, freeze it by renaming it <demoxx>
-
-
-%%  #1 parameter file name
-MAPparamsName='Normal';
-
-
-%% #2 probability (fast) or spikes (slow) representation
-AN_spikesOrProbability='spikes';
-% or
-AN_spikesOrProbability='probability';
-
-
-%% #3 pure tone, harmonic sequence or speech file input
-signalType= 'tones';
-duration=0.100;                 % seconds
-% duration=0.020;                 % seconds
-sampleRate= 64000;
-% toneFrequency= 250:250:8000;    % harmonic sequence (Hz)
-toneFrequency= 1000;            % or a pure tone (Hz8
-
-rampDuration=.005;              % seconds
-
-% or
-signalType= 'file';
-fileName='twister_44kHz';
-% fileName='new-da-44khz';
-
-% ? and mix with an optional second file?
-mixerFile=[];
-%or
-mixerFile='babble';
-leveldBSPL2=60;
-
-%% #4 rms level
-% signal details
-leveldBSPL= 70;                  % dB SPL
-
-
-%% #5 number of channels in the model
-%   21-channel model (log spacing)
-numChannels=21;
-lowestBF=250; 	highestBF= 8000;
-BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels));
-
-%   or specify your own channel BFs
-% BFlist=toneFrequency;
-
-
-%% #6 change model parameters
-paramChanges=[];
-
-% or
-% Parameter changes can be used to change one or more model parameters
-%  *after* the MAPparams file has been read
-% This example declares only one fiber type with a calcium clearance time
-% constant of 80e-6 s (HSR fiber) when the probability option is selected.
-% paramChanges={'AN_IHCsynapseParams.ANspeedUpFactor=5;', ...
-%     'IHCpreSynapseParams.tauCa=86e-6;'};
-% paramChanges={'DRNLParams.rateToAttenuationFactorProb = 0;'};
-
-%% delare showMap options
-showMapOptions=[];  % use defaults
-
-% or (example: show everything including an smoothed SACF output
-showMapOptions.showModelParameters=1;
-showMapOptions.showModelOutput=1;
-showMapOptions.printFiringRates=1;
-showMapOptions.showACF=0;
-showMapOptions.showEfferent=1;
-if strcmp(AN_spikesOrProbability, 'probability')
-    showMapOptions.surfProbability=1;
-end
-if strcmp(signalType, 'file')
-    showMapOptions.fileName=fileName;
-end
-
-%% Generate stimuli
-
-dbstop if error
-restorePath=path;
-addpath (['..' filesep 'MAP'],    ['..' filesep 'wavFileStore'])
-switch signalType
-    case 'tones'
-        inputSignal=createMultiTone(sampleRate, toneFrequency, ...
-            leveldBSPL, duration, rampDuration);
-        
-    case 'file'
-        %% file input simple or mixed
-        [inputSignal sampleRate]=wavread(fileName);
-            dt=1/sampleRate;
-        inputSignal=inputSignal(:,1);
-        targetRMS=20e-6*10^(leveldBSPL/20);
-        rms=(mean(inputSignal.^2))^0.5;
-        amp=targetRMS/rms;
-        inputSignal=inputSignal*amp;
-        silence= zeros(1,round(0.1/dt));
-        inputSignal= [silence inputSignal' silence];
-        
-        if ~isempty(mixerFile)
-            [inputSignal2 sampleRate2]=wavread(mixerFile);
-            if ~isequal(sampleRate,sampleRate2)
-                error...
-                    ('file and mixer file have different sample rates')
-            end
-            inputSignal2=inputSignal2(:,1);
-            [r c]=size(inputSignal);
-            inputSignal2=inputSignal2(1:r);
-            targetRMS=20e-6*10^(leveldBSPL2/20);
-            rms=(mean(inputSignal2.^2))^0.5;
-            amp=targetRMS/rms;
-            inputSignal2=inputSignal2*amp;
-            rampDuration=dt*length(inputSignal2)/3;
-            if rampDuration>0.5*duration, rampDuration=duration/2; end
-            rampTime=dt:dt:rampDuration;
-            time=dt*(1:length(inputSignal2));
-            ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
-                ones(1,length(time)-length(rampTime))];
-            inputSignal2=inputSignal2.*ramp';
-            ramp=fliplr(ramp);
-            inputSignal2=inputSignal2.*ramp';
-        end
-        
-        inputSignal=inputSignal+inputSignal2;
-end
-
-
-%% run the model
-tic
-
-fprintf('\n')
-disp(['Signal duration= ' num2str(length(inputSignal)/sampleRate)])
-disp([num2str(numChannels) ' channel model'])
-disp('Computing ...')
-
-restorePath=path;
-addpath (['..' filesep 'MAP'])
-
-MAP1_14(inputSignal, sampleRate, BFlist, ...
-    MAPparamsName, AN_spikesOrProbability, paramChanges);
-path(restorePath)
-toc
-
-% the model run is now complete. Now display the results
-showMAP(showMapOptions)
-for i=1:length(paramChanges)
-    disp(paramChanges{i})
-end
-
-toc
-path(restorePath)
-
-
-function inputSignal=createMultiTone(sampleRate, toneFrequency, ...
-    leveldBSPL, duration, rampDuration)
-% Create pure tone stimulus
-dt=1/sampleRate; % seconds
-time=dt: dt: duration;
-inputSignal=sum(sin(2*pi*toneFrequency'*time), 1);
-amp=10^(leveldBSPL/20)*28e-6;   % converts to Pascals (peak)
-inputSignal=amp*inputSignal;
-
-% apply ramps
-% catch rampTime error
-if rampDuration>0.5*duration, rampDuration=duration/2; end
-rampTime=dt:dt:rampDuration;
-ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
-    ones(1,length(time)-length(rampTime))];
-inputSignal=inputSignal.*ramp;
-ramp=fliplr(ramp);
-inputSignal=inputSignal.*ramp;
-
-% add 10 ms silence
-silence= zeros(1,round(0.03/dt));
-inputSignal= [silence inputSignal silence];
-
diff -r fafe69c43108 -r 45f28c49461e utilities/UTIL_PSTHmakerb.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/utilities/UTIL_PSTHmakerb.m	Mon Jun 13 18:21:05 2011 +0100
@@ -0,0 +1,37 @@
+function [PSTH ]=UTIL_PSTHmakerb(inputData, dt, PSTHbinWidth)
+% UTIL_PSTHmakerb averages mean values into bins.
+% No corrections are applied
+% usage:
+%	PSTH=UTIL_PSTHmaker(inputData, method)
+%
+% arguments
+%	inputData is a channel x time matrix
+%	PSTH is the reduced matrix, the sum of all elements in the bin
+%
+
+[numChannels numDataPoints]= size(inputData);
+
+% Multiple fibers is the same as repeat trials
+% Consolidate data into a histogram 
+dataPointsPerBin=round(PSTHbinWidth/dt);
+if dataPointsPerBin<=1;
+% 	Too few data points
+	PSTH=inputData;
+	return
+end
+
+numBins=floor(numDataPoints/dataPointsPerBin);
+PSTH=zeros(numChannels,numBins);
+
+% take care that signal length is an integer multiple of bin size
+%  by dropping the last unuseable values
+useableDataLength=numBins*dataPointsPerBin;
+inputData=inputData(:,1:useableDataLength);
+
+for ch=1:numChannels
+	% Convert each channel into a matrix where each column represents 
+	% the content of a single PSTH bin
+	PSTH2D=reshape (inputData(ch,:), dataPointsPerBin, numBins );
+	% and sum within each bin (across the rows
+	PSTH(ch,:)=mean (PSTH2D,1);
+end