Mercurial > hg > map
changeset 22:45f28c49461e master
removing duplicate changes
| author | Ray Meddis <rmeddis@essex.ac.uk> |
|---|---|
| date | Mon, 13 Jun 2011 18:21:05 +0100 |
| parents | c489ebada16e (diff) fafe69c43108 (current diff) |
| children | |
| files | MAP/MAP1_14.m old files/MAPparamsNormal.m old files/showMAP.m old files/test_MAP1_14.m parameterStore/MAPparamsNormal.m testPrograms/showMAP.m testPrograms/test_MAP1_14.m |
| diffstat | 7 files changed, 945 insertions(+), 822 deletions(-) [+] |
line wrap: on
line diff
--- /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 + +
--- /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)
--- /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]; +
--- 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 - -
--- 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)
--- 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]; -
--- /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