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 (current diff) fafe69c43108 (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 | 18 files changed, 986 insertions(+), 887 deletions(-) [+] |
line wrap: on
line diff
--- a/MAP/MAP1_14.m Mon Jun 13 18:13:29 2011 +0100 +++ b/MAP/MAP1_14.m Mon Jun 13 18:21:05 2011 +0100 @@ -535,34 +535,34 @@ % shorter segments after speed up. shorterSegmentEndPTR=reducedSegmentPTR+reducedSegmentLength-1; - iputPressureSegment=inputSignal... + inputPressureSegment=inputSignal... (:,segmentStartPTR:segmentStartPTR+segmentLength-1); % segment debugging plots % figure(98) - % plot(segmentTime,iputPressureSegment), title('signalSegment') + % plot(segmentTime,inputPressureSegment), title('signalSegment') % OME ---------------------- % OME Stage 1: external resonances. Add to inputSignal pressure wave - y=iputPressureSegment; + y=inputPressureSegment; for n=1:nOMEExtFilters % any number of resonances can be used [x OMEExtFilterBndry{n}] = ... filter(ExtFilter_b{n},ExtFilter_a{n},... - iputPressureSegment, OMEExtFilterBndry{n}); + inputPressureSegment, OMEExtFilterBndry{n}); x= x* OMEgainScalars(n); % This is a parallel resonance so add it y=y+x; end - iputPressureSegment=y; - OMEextEarPressure(segmentStartPTR:segmentEndPTR)= iputPressureSegment; + inputPressureSegment=y; + OMEextEarPressure(segmentStartPTR:segmentEndPTR)= inputPressureSegment; % OME stage 2: convert input pressure (velocity) to % tympanic membrane(TM) displacement using low pass filter [TMdisplacementSegment OME_TMdisplacementBndry] = ... - filter(TMdisp_b,TMdisp_a,iputPressureSegment, ... + filter(TMdisp_b,TMdisp_a,inputPressureSegment, ... OME_TMdisplacementBndry); % and save it TMoutput(segmentStartPTR:segmentEndPTR)= TMdisplacementSegment; @@ -610,6 +610,7 @@ else % no MOC available yet MOC=ones(1, segmentLength); end + plot(MOC) % current channel % first gammatone filter for order = 1 : GTnonlinOrder @@ -634,7 +635,8 @@ % second filter removes distortion products for order = 1 : GTnonlinOrder [ nonlinOutput GTnonlinBdry2{BFno,order}] = ... - filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), nonlinOutput, GTnonlinBdry2{BFno,order}); + filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ... + nonlinOutput, GTnonlinBdry2{BFno,order}); end % combine the two paths to give the DRNL displacement @@ -757,14 +759,18 @@ % Estimate efferent effects. ARattenuation (0 <> 1) % acoustic reflex - ARAttSeg=mean(ANrate(1:nBFs,:),1); %LSR channels are 1:nBF - % smooth - [ARAttSeg, ARboundaryProb] = ... - filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundaryProb); - ARAttSeg=ARAttSeg-ARrateThreshold; - ARAttSeg(ARAttSeg<0)=0; % prevent negative strengths - ARattenuation(segmentStartPTR:segmentEndPTR)=... - (1-ARrateToAttenuationFactorProb.* ARAttSeg); + [r c]=size(ANrate); + if r>nBFs % Only if LSR fibers are computed + ARAttSeg=mean(ANrate(1:nBFs,:),1); %LSR channels are 1:nBF + % smooth + [ARAttSeg, ARboundaryProb] = ... + filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundaryProb); + ARAttSeg=ARAttSeg-ARrateThreshold; + ARAttSeg(ARAttSeg<0)=0; % prevent negative strengths + ARattenuation(segmentStartPTR:segmentEndPTR)=... + (1-ARrateToAttenuationFactorProb.* ARAttSeg); + end + % plot(ARattenuation) % MOC attenuation % within-channel HSR response only @@ -781,6 +787,8 @@ end MOCattenuation(MOCattenuation<0)=0.001; + % plot(MOCattenuation) + case 'spikes' ANtimeCount=0;
--- a/multithreshold 1.46/expGUI_MT.m Mon Jun 13 18:13:29 2011 +0100 +++ b/multithreshold 1.46/expGUI_MT.m Mon Jun 13 18:21:05 2011 +0100 @@ -1292,8 +1292,6 @@ betweenRuns.observationCount= []; betweenRuns.timesOfFirstReversals= []; betweenRuns.bestThresholdTracks=''; -betweenRuns.bestThresholdMeanTracks=''; -betweenRuns.bestThresholdMedianTracks=''; betweenRuns.levelTracks=''; betweenRuns.responseTracks=''; betweenRuns.slopeKTracks= [];
--- a/multithreshold 1.46/nextStimulus.m Mon Jun 13 18:13:29 2011 +0100 +++ b/multithreshold 1.46/nextStimulus.m Mon Jun 13 18:21:05 2011 +0100 @@ -82,7 +82,7 @@ end set(handles.textMSG,'BackgroundColor','w', 'ForegroundColor', 'b') - + % Now the serious business of crafting and presenting the stimulus errormsg= stimulusMakeAndPlay (handles); @@ -314,6 +314,23 @@ targetLevel=-100; % no target end +% ----------------------------- calibration of sound output +calibrationCorrectiondB=stimulusParameters.calibrationdB; +if calibrationCorrectiondB<-50 + if maskerFrequency==targetFrequency + load 'calibrationFile' % calibrationFrequency calibrationAttenutation + idx=find(calibrationFrequency==targetFrequency); + if isempty(idx) + error('Calibration bty file; frequency not found') + else + calibrationCorrectiondB=calibrationAttenutation(idx) + end + else + error('calibration by file requested but masker frequency is not the same as target') + end +end + + % -------------------------------------- Checks on excessive signal level % clipping is relevant only for soundcard use (not modelling) @@ -329,7 +346,7 @@ switch experiment.ear case {'left', 'right', 'diotic',... 'dichotic', 'dioticLeft', 'dichoticRight'} - clippingLevel=91+stimulusParameters.calibrationdB; + clippingLevel=91+calibrationCorrectiondB; soundCardMinimum=clippingLevel-20*log10(2^24); otherwise clippingLevel=inf; @@ -366,7 +383,7 @@ withinRuns.forceThreshold=NaN; return end - + case 'targetLevel' upperLevel=stimulusParameters.WRVlimits(2); lowerLevel=stimulusParameters.WRVlimits(1); @@ -376,7 +393,7 @@ num2str(max(targetLevel, cueTargetLevel)) ... ') is too high ***']; withinRuns.forceThreshold=upperLevel; - withinRuns.forceThreshold=NaN; + withinRuns.forceThreshold=NaN; return end if max(targetLevel, cueTargetLevel)< lowerLevel @@ -384,7 +401,7 @@ num2str(max(targetLevel, cueTargetLevel)) ... ') is too low ***']; withinRuns.forceThreshold=lowerLevel; - withinRuns.forceThreshold=NaN; + withinRuns.forceThreshold=NaN; return end if max(targetLevel, cueTargetLevel)> clippingLevel @@ -392,10 +409,10 @@ num2str(max(targetLevel, cueTargetLevel)) ... ') is clipping ***']; withinRuns.forceThreshold=upperLevel; - withinRuns.forceThreshold=NaN; + withinRuns.forceThreshold=NaN; return end - end + end case 'maskerDuration' % this is odd! but harmless if max(maskerDuration, cueMaskerDuration)> ... @@ -493,7 +510,7 @@ backgroundType=stimulusParameters.backgroundType; switch stimulusParameters.backgroundType case {'noiseDich', 'pinkNoiseDich'} -% case 'Dich' + % case 'Dich' % dich means put the background in the ear opposite to the target backgroundType=backgroundType(1:end-4); switch targetEar @@ -503,7 +520,7 @@ backgroundEar=1; end otherwise -% case {'none','noise', 'pinkNoise', 'TEN','babble'} + % case {'none','noise', 'pinkNoise', 'TEN','babble'} backgroundEar=targetEar; end @@ -515,9 +532,9 @@ globalStimParams.dt=dt; stimulusParameters.dt=dt; % for use later -% calibration of sound output -correctiondB=stimulusParameters.calibrationdB; -globalStimParams.audioOutCorrection=10^(correctiondB/20); + + +globalStimParams.audioOutCorrection=10^(calibrationCorrectiondB/20); % the output will be reduced by this amount in stimulusCreate % i.e. audio=audio/globalStimParams.audioOutCorrection % A 91 dB level will yield a peak amp of 1 for calibration=0 @@ -680,7 +697,7 @@ 'backgrounds and maskers'... filesep 'ten.wav']; [tenNoise, FS]=wavread(fileName); - + tenNoise=resample(tenNoise, globalStimParams.FS, FS); stimComponents(backgroundEar,componentNo).type='file'; stimComponents(backgroundEar,componentNo).stimulus=tenNoise'; @@ -705,8 +722,8 @@ +stimulusParameters.maskerDuration; stimulusParameters.testTargetEnds=... stimulusParameters.testTargetBegins+withinRuns.variableValue; -% case 'SRT' -% set(handles.editdigitInput,'visible','off') + % case 'SRT' + % set(handles.editdigitInput,'visible','off') otherwise stimulusParameters.testTargetBegins=targetDelay; stimulusParameters.testTargetEnds=targetDelay+targetDuration;
--- a/multithreshold 1.46/printReport.m Mon Jun 13 18:13:29 2011 +0100 +++ b/multithreshold 1.46/printReport.m Mon Jun 13 18:21:05 2011 +0100 @@ -26,15 +26,16 @@ if experiment.saveData saveFileName=['savedData/' experiment.name '_' experiment.date '_' experiment.paradigm]; else - % overwrite existing file just in case + % save this data (just in case) saveFileName=['savedData/mostRecentResults']; end experiment.minElapsed=etime(clock, betweenRuns.timeNow)/60; - save(saveFileName, 'experiment', 'stimulusParameters', 'betweenRuns', 'withinRuns', 'statsModel', 'expGUIhandles') + save(saveFileName, 'experiment', 'stimulusParameters',... + 'betweenRuns', 'withinRuns', 'statsModel', 'expGUIhandles') disp(['data saved as: ' saveFileName]); else - % reprint request + % reprint request (i.e print out old data) printReportGuide.fileName=fileName; load(printReportGuide.fileName) saveFileName=printReportGuide.fileName; @@ -46,17 +47,18 @@ end fprintf('******** multithreshold ') + x=pwd; disp(['version ' x(end-3:end)]) fprintf('\nName:\t%s ', experiment.name) fprintf('\nparadigm:\t%s ', experiment.paradigm) fprintf('\nEar:\t%s ', experiment.ear) + method=experiment.threshEstMethod; if stimulusParameters.includeCue && ~strcmp(method(1:6),'2I2AFC') method=[method '/ withCue']; end fprintf('\nmethod:\t%s ', method) fprintf('\ndate:\t%s ', experiment.date) - fprintf('\n\n') if isempty(betweenRuns.thresholds) @@ -64,12 +66,11 @@ end % prepare results as matrices ready to print +% first sort the actual sequence into a more readable sequence [idx1, idx2, var1values, var2values]=... - sortVariables(betweenRuns.variableList1, betweenRuns.variableList2, betweenRuns.var1Sequence, betweenRuns.var2Sequence); + sortVariables(betweenRuns.variableList1, betweenRuns.variableList2,... + betweenRuns.var1Sequence, betweenRuns.var2Sequence); -% if strcmp(betweenRuns.variableName2, 'none') -% betweenRuns.variableName2=' '; -% end header1=betweenRuns.variableName1; header2=betweenRuns.variableName2; header1 = strrep(header1, 'none', ' '); % none is not a useful header @@ -77,7 +78,8 @@ headers=strvcat([header1 '/'], header2); disp('thresholds') -msg=printTabTable(sortTablesForPrinting(idx1,idx2,var1values,var2values, betweenRuns.thresholds), headers); +msg=printTabTable(sortTablesForPrinting(idx1,idx2,... + var1values,var2values, betweenRuns.thresholds), headers); addToMsg(msg,0) fprintf('\n') @@ -85,7 +87,6 @@ betweenRuns.levelTracks=betweenRuns.levelTracks(idx1); betweenRuns.responseTracks=betweenRuns.responseTracks(idx1); betweenRuns.bestThresholdTracks=betweenRuns.bestThresholdTracks(idx1); - betweenRuns.levelTracks=betweenRuns.levelTracks(idx2); betweenRuns.responseTracks=betweenRuns.responseTracks(idx2); betweenRuns.bestThresholdTracks=betweenRuns.bestThresholdTracks(idx2); @@ -136,21 +137,6 @@ header=strvcat(header, 'mean'); end - disp(' '); disp('threshold (mean) tracks starting from the first reversal') - for i=1:length(betweenRuns.bestThresholdTracks) - if printReportGuide.HorizontalTracks - end - printTabTable(betweenRuns.bestThresholdMeanTracks{i}); - header=strvcat(header, 'mean'); - end - disp(' '); disp('threshold tracks (median) starting from the first reversal') - for i=1:length(betweenRuns.bestThresholdMedianTracks) - if printReportGuide.HorizontalTracks - end - printTabTable(betweenRuns.bestThresholdTracks{i}); - header=strvcat(header, 'mean'); - end - end switch experiment.ear
--- a/multithreshold 1.46/subjGUI_MT.m Mon Jun 13 18:13:29 2011 +0100 +++ b/multithreshold 1.46/subjGUI_MT.m Mon Jun 13 18:21:05 2011 +0100 @@ -244,8 +244,6 @@ betweenRuns.catchTrials=[]; betweenRuns.timesOfFirstReversals=[]; betweenRuns.bestThresholdTracks=[]; -betweenRuns.bestThresholdMeanTracks=[]; -betweenRuns.bestThresholdMedianTracks=[]; betweenRuns.levelTracks=[]; betweenRuns.responseTracks=[]; betweenRuns.slopeKTracks=[]; @@ -1260,7 +1258,7 @@ switch experiment.threshEstMethod case {'MaxLikelihood', 'oneIntervalUpDown'} % last value in the list - threshold=withinRuns.meanEstTrack(end); +% threshold=withinRuns.meanEstTrack(end); threshold=withinRuns.thresholdEstimateTrack(end); stdev=NaN; @@ -1302,23 +1300,13 @@ % add variable length tracks to cell arrays if withinRuns.beginningOfPhase2>0 betweenRuns.bestThresholdTracks{length(betweenRuns.thresholds)}=... - withinRuns.thresholdEstimateTrack; - betweenRuns.bestThresholdMeanTracks... - {length(betweenRuns.thresholds_mean)}=... - withinRuns.thresholdEstimateTrack; - betweenRuns.bestThresholdMedianTracks... - {length(betweenRuns.thresholds_median)}=... - withinRuns.thresholdEstimateTrack; - + withinRuns.thresholdEstimateTrack; betweenRuns.levelTracks{length(betweenRuns.thresholds)}=... withinRuns.levelList(withinRuns.beginningOfPhase2:end); betweenRuns.responseTracks{length(betweenRuns.thresholds)}=... withinRuns.responseList(withinRuns.beginningOfPhase2:end); else betweenRuns.bestThresholdTracks{length(betweenRuns.thresholds)}=[]; - betweenRuns.bestThresholdMeanTracks{length(betweenRuns.thresholds)}=[]; - betweenRuns.bestThresholdMedianTracks{length(betweenRuns.thresholds)}=... - []; betweenRuns.levelTracks{length(betweenRuns.thresholds)}=[]; betweenRuns.responseTracks{length(betweenRuns.thresholds)}=[]; end
--- a/multithreshold 1.46/testOME.m Mon Jun 13 18:13:29 2011 +0100 +++ b/multithreshold 1.46/testOME.m Mon Jun 13 18:21:05 2011 +0100 @@ -23,9 +23,11 @@ set(2,'position',[5 349 268 327]) semilogx(HuberFrequencies, 20*log10(HuberDisplacementAt80dBSPL/1e-10),... 'ko', 'MarkerFaceColor','k', 'Marker','o', 'markerSize',6) -% Generate test stimulus ................................................................. +hold on -%% independent test using discrete frequencies +%% Generate test stimulus ................................................................. + +% independent test using discrete frequencies peakResponses=[]; peakTMpressure=[]; frequencies=[200 400 HuberFrequencies 10000]; @@ -39,7 +41,7 @@ paramChanges{1}='OMEParams.rateToAttenuationFactorProb=0;'; paramChanges{2}='DRNLParams.rateToAttenuationFactorProb = 0;'; - global OMEoutput OMEextEarPressure TMoutput ARAttenuation + global OMEoutput OMEextEarPressure TMoutput ARattenuation % BF is irrelevant MAP1_14(inputSignal, sampleRate, -1, ... MAPparamsName, AN_spikesOrProbability, paramChanges); @@ -48,7 +50,7 @@ peakResponses=[peakResponses peakDisplacement]; peakTMpressure=[peakTMpressure max(OMEextEarPressure)]; - disp([' greatest AR attenuation: ' num2str(min(ARAttenuation))]) + disp([' AR attenuation (dB): ' num2str(20*log10(min(ARattenuation)))]) end %% Report @@ -59,6 +61,7 @@ % stapes peak displacement figure(2), subplot(2,1,1), hold on semilogx(frequencies, 20*log10(peakResponses/1e-10), 'r', 'linewidth',4) +set(gca,'xScale','log') % ylim([1e-11 1e-8]) xlim([100 10000]), ylim([0 30]) grid on
--- /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 18:13:29 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,340 +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.003; % * 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=3e2; % 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 = .007; % strength of MOC -DRNLParams.rateToAttenuationFactorProb = .002; % strength of MOC -% DRNLParams.rateToAttenuationFactorProb = .0; % strength of MOC -DRNLParams.MOCtau =.03; % 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=50e-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 18:13:29 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,273 +0,0 @@ -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)))]) - 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 - -%% 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 - -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) - 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 18:13:29 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,197 +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= 100000; -% toneFrequency= 250:250:8000; % harmonic sequence (Hz) -toneFrequency= 2000; % or a pure tone (Hz8 - -rampDuration=.005; % seconds - -% % or -% signalType= 'file'; -% fileName='twister_44kHz'; -% % fileName='new-da-44khz'; - -% mix with an optional second file? -mixerFile=[]; -%or -% mixerFile='babble'; -% leveldBSPL2=30; - -%% #4 rms level -% signal details -leveldBSPL= 60; % 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={'AN_IHCsynapseParams.ANspeedUpFactor=5;', ... -% '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); - inputSignal=inputSignal(:,1); - targetRMS=20e-6*10^(leveldBSPL/20); - rms=(mean(inputSignal.^2))^0.5; - amp=targetRMS/rms; - inputSignal=inputSignal*amp; - silence=zeros(10000,1); - inputSignal=[silence; inputSignal]; - if ~isempty(mixerFile) - [inputSignal2 sampleRate]=wavread(mixerFile); - 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; - inputSignal=inputSignal+inputSignal2; - end -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/utilities/UTIL_makePSTH.m Mon Jun 13 18:13:29 2011 +0100 +++ b/utilities/UTIL_makePSTH.m Mon Jun 13 18:21:05 2011 +0100 @@ -1,4 +1,5 @@ -function PSTH =UTIL_makePSTH(inputData, dt, PSTHbinWidth) +function [PSTH dataPointsPerBin]= ... + UTIL_makePSTH(inputData, dt, PSTHbinWidth) % UTIL_PSTHmaker accumulates values into bins. % No corrections are applied % usage: