Mercurial > hg > map
changeset 29:b51bf546ca3f
physiologyProb
line wrap: on
line diff
--- a/multithreshold 1.46/expGUI_MT.m Fri Jul 01 12:59:47 2011 +0100 +++ b/multithreshold 1.46/expGUI_MT.m Fri Jul 08 13:48:27 2011 +0100 @@ -1,9 +1,9 @@ % expGUI_MT = 'experimenter GUI for multiThreshold % allows the experimenter to design experiments. % The *running* of experiments is left to subjGUI.m -% +% % There are three kinds of experiments known as 'earOptions': -% 1. Measurements using real listeners +% 1. Measurements using real listeners % 'left', 'right', 'diotic', 'dichoticLeft', 'dichoticRight' % 2. Measurements using the MAP model as the subject % 'MAPmodelListen', 'MAPmodelMultiCh', 'MAPmodelSingleCh' @@ -88,7 +88,7 @@ setLocationOfGUIs(handles) function setLocationOfGUIs(handles) -global checkForPreviousGUI % holds screen positioning across repeated calls +global checkForPreviousGUI % holds screen positioning across repeated calls scrnsize=get(0,'screensize'); checkForPreviousGUI=[]; % if isstruct(checkForPreviousGUI)... @@ -140,6 +140,7 @@ % addpath (['..' filesep 'modules'], ['..' filesep 'utilities'], ... % ['..' filesep 'parameterStore'], ['..' filesep 'wavFileStore'],... % ['..' filesep 'testPrograms']) +addpath (['..' filesep 'testPrograms']) % specify all variables that need to be set on the GUI variableNames={'stimulusDelay','maskerDuration','maskerLevel',... @@ -314,7 +315,7 @@ paradigm=paradigmNames{chosenOption}; experiment.paradigm=paradigm; - %Paradigm: read in all relevant parameters +%Paradigm: read in all relevant parameters % a file must exist with this name 'paradigm_<paradigm>' % 'handles' are only occasionally used addpath ('paradigms') @@ -347,7 +348,7 @@ eval (cmd); % establish popup menus on the basis of the paradigm file -set(handles.popupmenuRandomize,'value', betweenRuns.randomizeSequence) +set(handles.popupmenuRandomize,'value', betweenRuns.randomizeSequence) set(handles.popupmenuPhase,'string', stimulusParameters.maskerPhase) if stimulusParameters.includeCue set(handles.popupmenuCueNoCue,'value', 1) @@ -383,14 +384,14 @@ % ------------------------------------------------------ aResetPopupMenus function aResetPopupMenus(handles) -global stimulusParameters betweenRuns variableNames +global stimulusParameters betweenRuns variableNames global targetTypes maskerTypes experiment backgroundTypes switch experiment.threshEstMethod case {'MaxLikelihood','oneIntervalUpDown'} set(handles.editstopCriteriaBox, 'string', ... num2str(experiment.singleIntervalMaxTrials)) - + case {'2I2AFC++','2I2AFC+++'} set(handles.editstopCriteriaBox, 'string', ... num2str(experiment.stopCriteria2IFC)) @@ -410,7 +411,7 @@ for i=1:length(variableNames) if strcmp(variableNames{i},betweenRuns.variableName1) variableParameter1ID=i; - end + end if strcmp(variableNames{i},betweenRuns.variableName2) variableParameter2ID=i; end @@ -453,7 +454,7 @@ set(handles.editStatsModel, 'visible', 'on') set(handles.textStatsModel, 'visible', 'on') set(handles.pushbuttonStop, 'visible', 'on') -showModelPushButtons(handles, 0) + showModelPushButtons(handles, 0) set(handles.editCatchTrialRate, 'visible', 'off') set(handles.textCatchTrials, 'visible', 'off') set(handles.editcalibrationdB, 'visible', 'off') @@ -462,25 +463,25 @@ set(handles.textCue, 'visible', 'off') set(handles.editMusicLevel,'visible', 'off') set(handles.textMusicLevel,'visible', 'off') - + case {'MAPmodel', 'MAPmodelListen', 'MAPmodelMultiCh', 'MAPmodelSingleCh'} set(handles.popupmenuCueNoCue, 'visible', 'off') set(handles.editStatsModel, 'visible', 'off') set(handles.textStatsModel, 'visible', 'off') set(handles.pushbuttonStop, 'visible', 'on') -showModelPushButtons(handles, 1) + showModelPushButtons(handles, 1) set(handles.editcalibrationdB, 'visible', 'off') set(handles.textcalibration, 'visible', 'off') set(handles.textCue, 'visible', 'off') set(handles.editMusicLevel,'visible', 'off') set(handles.textMusicLevel,'visible', 'off') - + otherwise % i.e. using real subjects (left, right, diotic, dichotic) set(handles.editStatsModel, 'visible', 'off') set(handles.textStatsModel, 'visible', 'off') set(handles.pushbuttonStop, 'visible', 'off') -showModelPushButtons(handles, 0) + showModelPushButtons(handles, 0) set(handles.editCatchTrialRate, 'visible', 'on') set(handles.textCatchTrials, 'visible', 'on') set(handles.editcalibrationdB, 'visible', 'on') @@ -497,7 +498,7 @@ set(handles.textCue,'visible', 'on') set(handles.editstopCriteriaBox, 'string', ... num2str(experiment.singleIntervalMaxTrials)) - + if stimulusParameters.includeCue==0 set(handles.editcueTestDifference,'visible', 'off') set(handles.textcueTestDifference,'visible', 'off') @@ -505,7 +506,7 @@ set(handles.editcueTestDifference,'visible', 'on') set(handles.textcueTestDifference,'visible', 'on') end - + case {'2I2AFC++','2I2AFC+++'} set(handles.editCatchTrialRate, 'visible', 'off') set(handles.textCatchTrials, 'visible', 'off') @@ -609,7 +610,7 @@ % ------------------------------------------------ pushbuttonRun_Callback function pushbuttonRun_Callback(hObject, eventdata, handles) -global checkForPreviousGUI % holds screen positioning across repeated calls +global checkForPreviousGUI % holds screen positioning across repeated calls global experiment betweenRuns paradigmNames errormsg checkForPreviousGUI.GUIposition=get(handles.figure1,'position'); experiment.singleShot=0; @@ -621,8 +622,8 @@ set(handles.edittargetDuration,'string', num2str(0.25)) set(handles.editstopCriteriaBox,'string','10') % nTrials run (handles) - - if ~isempty(errormsg) + + if strcmp(errormsg,'manually stopped') disp(errormsg) optionNo=strmatch('profile',paradigmNames); set(handles.popupmenuParadigm,'value',optionNo); @@ -630,14 +631,14 @@ aParadigmSelection(handles) return end - + global resultsTable longTone=resultsTable(2:end,2:end); - + set(handles.edittargetDuration,'string', num2str(0.016)) set(handles.editstopCriteriaBox,'string','10') % nTrials run (handles) - if ~isempty(errormsg) + if strcmp(errormsg,'manually stopped') disp(errormsg) optionNo=strmatch('profile',paradigmNames); set(handles.popupmenuParadigm,'value',optionNo); @@ -646,7 +647,7 @@ aParadigmSelection(handles) return end - + shortTone=resultsTable(2:end,2:end); % use these threshold for TMC @@ -656,10 +657,10 @@ aParadigmSelection(handles) set(handles.edittargetLevel,'string', thresholds16ms+10); set(handles.editstopCriteriaBox,'string','10') % nTrials - pause(.1) + pause(.1) run (handles) - - if ~isempty(errormsg) + + if strcmp(errormsg,'manually stopped') disp(errormsg) optionNo=strmatch('profile',paradigmNames); set(handles.popupmenuParadigm,'value',optionNo); @@ -672,17 +673,29 @@ TMC=resultsTable(2:end,2:end); gaps=resultsTable(2:end,1); BFs=resultsTable(1, 2:end); - + % use these threshold for IFMC optionNo=strmatch('IFMC',paradigmNames); set(handles.popupmenuParadigm,'value',optionNo); aParadigmSelection(handles) set(handles.edittargetLevel,'string', thresholds16ms+10); set(handles.editstopCriteriaBox,'string','10') % nTrials - pause(.1) + pause(.1) run (handles) - if ~isempty(errormsg) + + IFMCs=resultsTable(2:end,2:end); + offBFs=resultsTable(2:end,1); + + % reset original paradigm + optionNo=strmatch('profile',paradigmNames); + set(handles.popupmenuParadigm,'value',optionNo); + aParadigmSelection(handles) + + save profile longTone shortTone gaps BFs TMC offBFs IFMCs + plotProfile(longTone,shortTone,gaps,BFs,TMC,offBFs,IFMCs) + + if strcmp(errormsg,'manually stopped') disp(errormsg) optionNo=strmatch('profile',paradigmNames); set(handles.popupmenuParadigm,'value',optionNo); @@ -690,19 +703,7 @@ experiment.stop=-0; aParadigmSelection(handles) return - end - - IFMCs=resultsTable(2:end,2:end); - offBFs=resultsTable(2:end,1); - - % reset original paradigm - optionNo=strmatch('profile',paradigmNames); - set(handles.popupmenuParadigm,'value',optionNo); - aParadigmSelection(handles) - - save profile longTone shortTone gaps BFs TMC offBFs IFMCs - plotProfile(longTone,shortTone,gaps,BFs,TMC,offBFs,IFMCs) - + end otherwise run (handles) experiment.stop=0; @@ -713,7 +714,7 @@ tic expGUIhandles=handles; set(handles.pushbuttonStop, 'backgroundColor', [.941 .941 .941]) -set(handles.editparamChanges,'visible','off') +% set(handles.editparamChanges,'visible','off') % message box white (removes any previous error message) set(handles.textMSG,... @@ -742,7 +743,7 @@ % --- Executes on button press in pushbuttonSingleShot. function pushbuttonSingleShot_Callback(hObject, eventdata, handles) -global experiment +global experiment experiment.singleShot=1; % special test for spontaneous activity @@ -764,7 +765,7 @@ % ------------------------------------------aReadAndCheckParameterBoxes function errorMsg=aReadAndCheckParameterBoxes(handles) global experiment stimulusParameters betweenRuns statsModel -global variableNames LevittControl +global variableNames LevittControl paramChanges % When the program sets the parameters all should be well % But when the user changes them... @@ -780,8 +781,8 @@ switch experiment.ear case { 'MAPmodel', 'MAPmodelMultiCh', ... - 'MAPmodelSingleCh', 'MAPmodelListen'} - % MAPmodel writes forced parameter settings to the screen + 'MAPmodelSingleCh', 'MAPmodelListen'} + % MAPmodel writes forced parameter settings to the screen % so that they can be read from there set(handles.popupmenuRandomize,'value',2) % fixed sequence set(handles.editstimulusDelay,'string','0.01') % no stimulus delay @@ -840,13 +841,13 @@ end % calibration -% this is used to *reduce* the output signal from what it otherwise +% this is used to *reduce* the output signal from what it otherwise % would be % signal values are between 1 - 2^23 % these are interpreted as microPascals between -29 dB and 128 dB SPL -% calibrationdB adjusts these values to compensate for equipment +% calibrationdB adjusts these values to compensate for equipment % characteristics -% this will change the range. e.g. a 7 dB calibration will yield +% this will change the range. e.g. a 7 dB calibration will yield % a range of -36 to 121 dB SPL % Calibration is not used when modelling. Values are treated as dB SPL stimulusParameters.calibrationdB=... @@ -946,7 +947,7 @@ % start value for step until reduced LevittControl.startLevelStep= stimulusParameters.WRVsteps(1); % reduced step size - LevittControl.steadyLevittStep= stimulusParameters.WRVsteps(2); + LevittControl.steadyLevittStep= stimulusParameters.WRVsteps(2); LevittControl.TurnsToSmallSteps= 2; LevittControl.useLastNturns= 2*experiment.peaksUsed; LevittControl.minReversals= ... @@ -1006,7 +1007,7 @@ % Check WRVstartValues for length and compatibility with randomization % only one start value supplied so all start values are the same -if length(stimulusParameters.WRVstartValues)==1 +if length(stimulusParameters.WRVstartValues)==1 stimulusParameters.WRVstartValues= ... repmat(stimulusParameters.WRVstartValues, 1, ... length(betweenRuns.variableList1)... @@ -1051,7 +1052,7 @@ stimulusParameters.maskerDuration) addToMsg(... 'Warning: masker and target duration not the same.',1,1) - end + end if ~isequal(stimulusParameters.maskerLevel, ... stimulusParameters.targetLevel) addToMsg(['Warning: masker and target level different.'... @@ -1059,6 +1060,10 @@ end end +% identify model parameter changes if any +paramChanges=get(handles.editparamChanges,'string'); +eval(paramChanges); + % -------------------------------------------- aSetSampleRate function aSetSampleRate(sampleRate, handles) global stimulusParameters @@ -1106,24 +1111,24 @@ case {'statsModelLogistic', 'statsModelRareEvent'} set(handles.editStatsModel, 'visible', 'off') set(handles.textStatsModel, 'visible', 'off') - + % default psychometric bin size and logistic slopes set(handles.popupmenuRandomize,'value',2) % fixed sequence set(handles.editName,'string', 'statsModel') % experiment.headphonesUsed=0; - + case {'MAPmodelListen', 'MAPmodelMultiCh', 'MAPmodelSingleCh'} stimulusParameters.includeCue=0; % no cue set(handles.popupmenuCueNoCue,'value', 2) - + set(handles.editCatchTrialRate,'string','0 0 2 ');%no catch trials set(handles.editName,'string', 'Normal') % force name experiment.name=get(handles.editName,'string'); % read name back - set(handles.editcalibrationdB,'string','0') - + set(handles.editcalibrationdB,'string','0') + set(handles.popupmenuRandomize,'value',2) % fixed sequence set(handles.editstimulusDelay,'string','0') - + set(handles.editSaveData,'string', '0') set(handles.editSubjectFont,'string', '10'); experiment.MacGThreshold=0; % num MacG spikes to exceed threshold @@ -1178,7 +1183,7 @@ case 'training' experiment.possLogSlopes=0.5; end - + case 'oneIntervalUpDown' experiment.functionEstMethod='logisticLS'; set(handles.textstopCriteria,'string', 'stop criteria \ maxTrials') @@ -1189,7 +1194,7 @@ otherwise experiment.allowCatchTrials= 1; end - + case {'2I2AFC++', '2I2AFC+++'} LevittControl.rule='++'; % e.g. '++' or '+++' experiment.singleIntervalMaxTrials=experiment.stopCriteria2IFC; @@ -1209,7 +1214,7 @@ % NB responsibility for this is now transferred to the paradigm file switch experiment.threshEstMethod % only one value required for level change - case {'2I2AFC++', '2A2AIFC+++'} + case {'2I2AFC++', '2A2AIFC+++'} stimulusParameters.subjectText=... 'did the tone occur in window 1 or 2?'; case {'MaxLikelihood', 'oneIntervalUpDown'}; @@ -1399,41 +1404,49 @@ drawnow function pushbuttonOME_Callback(hObject, eventdata, handles) -global experiment +global experiment paramChanges aReadAndCheckParameterBoxes(handles); -testOME(experiment.name); +testOME(experiment.name, paramChanges); function pushbuttonBM_Callback(hObject, eventdata, handles) -global stimulusParameters experiment +global stimulusParameters experiment paramChanges aReadAndCheckParameterBoxes(handles); relativeFrequencies=[0.25 .5 .75 1 1.25 1.5 2]; - +AN_spikesOrProbability='probability'; testBM(stimulusParameters.targetFrequency, ... - experiment.name,relativeFrequencies); + experiment.name,relativeFrequencies, AN_spikesOrProbability, ... + paramChanges); function pushbuttonAN_Callback(hObject, eventdata, handles) -global stimulusParameters +global stimulusParameters experiment paramChanges aReadAndCheckParameterBoxes(handles); % now carry out tests showPSTHs=0; targetFrequency=stimulusParameters.targetFrequency(1); BFlist=targetFrequency; -testAN(targetFrequency,BFlist); +testAN(targetFrequency,BFlist,-10:10:90,experiment.name, paramChanges); function pushbuttonPhLk_Callback(hObject, eventdata, handles) +global experiment aReadAndCheckParameterBoxes(handles); -testPhaseLocking +testPhaseLocking(experiment.name) function pushbuttonSYN_Callback(hObject, eventdata, handles) +global stimulusParameters experiment paramChanges aReadAndCheckParameterBoxes(handles); -testSynapse +% now carry out tests +showPSTHs=0; +targetFrequency=stimulusParameters.targetFrequency(1); +BFlist=targetFrequency; +testSynapse(BFlist,experiment.name, paramChanges) function pushbuttonFM_Callback(hObject, eventdata, handles) -global stimulusParameters experiment +global stimulusParameters experiment paramChanges aReadAndCheckParameterBoxes(handles); showPSTHs=1; -testFM(stimulusParameters.targetFrequency(1),experiment.name, showPSTHs) +testFM(stimulusParameters.targetFrequency(1),experiment.name, ... + showPSTHs, paramChanges) function popupmenuPhase_Callback(hObject, eventdata, handles) global stimulusParameters @@ -1447,10 +1460,9 @@ function pushbuttonParams_Callback(hObject, eventdata, handles) global experiment stimulusParameters aReadAndCheckParameterBoxes(handles); -% print model parameters using the 'name' box (e.g. CTa -> MAPparamsCTa) showParams=1; BFlist=-1; paramChanges=get(handles.editparamChanges,'string'); -eval(paramChanges) +eval(paramChanges); paramFunctionName=['method=MAPparams' experiment.name ... '(BFlist, stimulusParameters.sampleRate, showParams,paramChanges);']; @@ -1459,10 +1471,10 @@ % --- Executes on button press in pushbuttonRP. function pushbuttonRP_Callback(hObject, eventdata, handles) -global experiment stimulusParameters +global experiment stimulusParameters paramChanges aReadAndCheckParameterBoxes(handles); % now carry out test -testRP(stimulusParameters.targetFrequency,experiment.name) +testRP(stimulusParameters.targetFrequency,experiment.name, paramChanges) % function handles % ?? @@ -1755,7 +1767,7 @@ end - + function editparamChanges_Callback(hObject, eventdata, handles)
--- a/multithreshold 1.46/nextStimulus.m Fri Jul 01 12:59:47 2011 +0100 +++ b/multithreshold 1.46/nextStimulus.m Fri Jul 08 13:48:27 2011 +0100 @@ -7,13 +7,13 @@ errormsg=''; % interrupt by 'stop' button -if experiment.stop - disp('******** experiment manually stopped *****************') - experiment.status= 'waitingForStart'; - errormsg='manually stopped'; - addToMsg(errormsg,1) - return -end +% if experiment.stop +% disp('******** experiment manually stopped *****************') +% experiment.status= 'waitingForStart'; +% errormsg='manually stopped'; +% addToMsg(errormsg,1) +% return +% end % -----------------------------------------choose catch trials at random % catch trials are for subject threshold measurements only
--- a/multithreshold 1.46/old files/MAPmodel.m Fri Jul 01 12:59:47 2011 +0100 +++ b/multithreshold 1.46/old files/MAPmodel.m Fri Jul 08 13:48:27 2011 +0100 @@ -4,7 +4,7 @@ global outerMiddleEarParams DRNLParams AN_IHCsynapseParams savePath=path; -addpath('..\MAP') +addpath(['..' filesep 'MAP'], ['..' filesep 'utilities']) modelResponse=[]; MacGregorResponse=[]; @@ -25,12 +25,12 @@ MAPparamsName, AN_spikesOrProbability); if showPlotsAndDetails - options.showModelParameters=0; + options.printModelParameters=0; options.showModelOutput=1; options.printFiringRates=1; options.showACF=0; options.showEfferent=1; - UTIL_showMAP(options, paramChanges) + UTIL_showMAP(options) end % No response, probably caused by hitting 'stop' button
--- a/multithreshold 1.46/paradigms/paradigm_IFMC.m Fri Jul 01 12:59:47 2011 +0100 +++ b/multithreshold 1.46/paradigms/paradigm_IFMC.m Fri Jul 08 13:48:27 2011 +0100 @@ -4,7 +4,7 @@ paradigm_training(handles) % default stimulusParameters.WRVname='maskerLevel'; -stimulusParameters.WRVstartValues=-10; +stimulusParameters.WRVstartValues=50; stimulusParameters.WRVsteps= [-10 -2]; stimulusParameters.WRVlimits=[-30 110];
--- a/multithreshold 1.46/paradigms/paradigm_TMC.m Fri Jul 01 12:59:47 2011 +0100 +++ b/multithreshold 1.46/paradigms/paradigm_TMC.m Fri Jul 08 13:48:27 2011 +0100 @@ -4,7 +4,7 @@ paradigm_training(handles) % default stimulusParameters.WRVname='maskerLevel'; -stimulusParameters.WRVstartValues=-10; +stimulusParameters.WRVstartValues=50; stimulusParameters.WRVsteps= [-10 -4]; stimulusParameters.WRVlimits=[-30 110];
--- a/multithreshold 1.46/paradigms/paradigm_profile.m Fri Jul 01 12:59:47 2011 +0100 +++ b/multithreshold 1.46/paradigms/paradigm_profile.m Fri Jul 08 13:48:27 2011 +0100 @@ -5,7 +5,6 @@ betweenRuns.variableName1='targetFrequency'; betweenRuns.variableList1=[250 500 1000 2000 4000 8000 ]; -betweenRuns.variableList1=str2num(get(handles.edittargetFrequency,'string')); betweenRuns.variableName2='targetDuration'; betweenRuns.variableList2=0.016; betweenRuns.randomizeSequence=1; % 'random sequence' @@ -15,7 +14,7 @@ stimulusParameters.targetType='tone'; stimulusParameters.targetPhase='sin'; -stimulusParameters.targetFrequency=1000; +stimulusParameters.targetFrequency=[250 500 1000 2000 4000 8000 ]; stimulusParameters.targetDuration=0.016; stimulusParameters.targetLevel=stimulusParameters.WRVstartValues(1);
--- a/multithreshold 1.46/plotProfile.m Fri Jul 01 12:59:47 2011 +0100 +++ b/multithreshold 1.46/plotProfile.m Fri Jul 08 13:48:27 2011 +0100 @@ -5,63 +5,118 @@ load profile end -normLongTone=[ 11.4 1.55 -13.5 -6.35 -6.4 7.45]; -normShortTone=[ 23.85 18.9 9.85 10.6 9.55 21.9]; +% comparison data (e.g. participants) +% rows are BFs -normGaps=0.01:0.01:0.09; -normTMC= [ -28.5 35.0 49.3 70.1 80.5 85.5 NaN NaN NaN; -31.1 44.3 48.4 59.5 56.4 76.7 70.2 82.4 76.3; -33.4 38.4 48.8 55.8 64.5 78.7 84.2 88.3 90.3; -25.4 37.0 49.2 49.7 58.2 69.6 87.7 95.8 93.0; -18.2 23.5 27.4 41.5 64.3 82.1 86.7 91.2 NaN; -32.5 35.8 43.5 52.1 69.1 78.6 86.6 86.0 NaN; - ]; -normTMC=normTMC'; +% -------------------------------------------JSan +compareBFs=[250 500 1000 2000 3000 +]; +compareLongTone= [67.46485 56.95655 65.01985 61.46655 73.33265 +]; +compareShortTone=[ 72.3185 63.2818 69.0373 65.2853 76 +]; -normIFMC=[ -50 42 34 35 34 33 37; -58 51 38 33 28 41 49; -57 41 27 20 28 37 66; -61 49 27 20 34 68 79; -67 45 27 22 46 74 87; -62 62 43 22 47 56 83; +compareGaps=[0.02 0.05 0.08]; +compareTMC= [ + 84 69 77 75 93 +88 73 81 79 97 +95 79 85 83 98 ]; -normIFMC=normIFMC'; + +compareMaskerFreqs=[0.7 0.9 1 1.1 1.3 ]; +compareIFMCs=[83.1698 77.3165 79.8474 82.9074 82.3294 +80.9667 73.6653 80.9446 80.7005 79.0022 +82.0135 71.2284 78.7345 74.3342 84.126 +79.3348 70.3347 78.5769 79.7274 91.9849 +79.2308 76.3164 83.6881 86.2105 NaN + ]; + +% % -------------------------------------------JE +% BFs=[250 500 1000 2000 4000 8000]; +% compareLongTone= [32 30 31 40 54 NaN]; +% compareShortTone=[ 49 50 47 56 63 NaN]; +% +% compareGaps=0.01:0.01:0.09; +% compareTMC= [ +% 69 83 82 NaN NaN NaN NaN NaN NaN +% 61 68 79 88 93 NaN NaN NaN NaN +% 63 69 79 84 92 NaN NaN NaN NaN +% 67 71 75 80 82 84 88 93 NaN +% 82 82 86 86 NaN 83 88 90 75 +% NaN NaN NaN NaN NaN NaN NaN NaN NaN +% ]; +% compareTMC=compareTMC'; +% +% compareMaskerFreqs=[0.5 0.7 0.9 1 1.1 1.3 1.6]; +% compareIFMCs=[ +% 64 60 58 58 57 60 64 +% 65 63 58 55 54 59 69 +% 68 64 60 59 62 73 79 +% 76 75 71 67 68 71 77 +% 79 71 68 69 73 75 77 +% 76 73 75 75 76 80 NaN +% ]; +% compareIFMCs=compareIFMCs'; + +% -------------------------------------------CMR +% CMR +% BFs=[250 500 1000 2000 4000 8000]; +% compareLongTone=[ 11.4 1.55 -13.5 -6.35 -6.4 7.45]; +% compareShortTone=[ 23.85 18.9 9.85 10.6 9.55 21.9]; +% +% compareGaps=0.01:0.01:0.09; +% compareTMC= [ +% 28.5 35.0 49.3 70.1 80.5 85.5 NaN NaN NaN; +% 31.1 44.3 48.4 59.5 56.4 76.7 70.2 82.4 76.3; +% 33.4 38.4 48.8 55.8 64.5 78.7 84.2 88.3 90.3; +% 25.4 37.0 49.2 49.7 58.2 69.6 87.7 95.8 93.0; +% 18.2 23.5 27.4 41.5 64.3 82.1 86.7 91.2 NaN; +% 32.5 35.8 43.5 52.1 69.1 78.6 86.6 86.0 NaN; +% ]; +% compareTMC=compareTMC'; +% +% compareMaskerFreqs=[0.5 0.7 0.9 1 1.1 1.3 1.6]; +% compareIFMCs=[ +% 50 42 34 35 34 33 37; +% 58 51 38 33 28 41 49; +% 57 41 27 20 28 37 66; +% 61 49 27 20 34 68 79; +% 67 45 27 22 46 74 87; +% 62 62 43 22 47 56 83; +% ]; +% compareIFMCs=compareIFMCs'; % absolute thresholds figure(90), clf subplot(2,1,2) semilogx(BFs,longTone,'ko-','lineWidth',2); hold on semilogx(BFs,shortTone,'bo-','lineWidth',2); hold on -semilogx(BFs,normLongTone,'ko:'); hold on -semilogx(BFs,normShortTone,'bo:'); hold on +semilogx(compareBFs,compareLongTone,'ko:'); hold on +semilogx(compareBFs,compareShortTone,'bo:'); hold on ylim([0 100]) % TMC for BFno=1:length(BFs) subplot(2,6,BFno) plot(gaps,TMC(:,BFno)-longTone(BFno),'r','lineWidth',3), hold on - plot(normGaps,normTMC(:,BFno)-longTone(BFno),'k:') - ylim([-10 90]) + plot(gaps,TMC(:,BFno),'b','lineWidth',3), hold on + ylim([-10 110]) xlim([0.01 0.1]) - grid on + grid on if BFno==1 ylabel('masker dB SL') xlabel('gap') - text(0.02,80,' TMC','backgroundColor','w') +% text(0.02,80,' TMC','backgroundColor','w') end - title([num2str(BFs(BFno)) ' Hz']) - set(gca,'XTick',[ 0.1],'xTickLabel', { '0.1'}) + title([num2str(BFs(BFno)) ' Hz']) + set(gca,'XTick',[ 0.1],'xTickLabel', { '0.1'}) end % IFMCs for BFno=1:length(BFs) - BF=BFs(BFno); freq=offBFs'*BFs(BFno); subplot(2,1,2) semilogx(freq,IFMCs(:,BFno),'r','lineWidth',3), hold on - semilogx(freq,normIFMC(:,BFno),'k:') ylim([0 100]) xlim([100 12000]) grid on @@ -69,3 +124,29 @@ xlabel('frequency (Hz)') ylabel('masker dB / probe dB') set(gca,'XTick',BFs) + +for BFno=1:length(compareBFs) + subplot(2,6,BFno) + plot(compareGaps,compareTMC(:,BFno)-longTone(BFno),'k:') + plot(compareGaps,compareTMC(:,BFno),'k:') + ylim([-10 110]) + xlim([0.01 0.1]) + grid on + if BFno==1 + ylabel('masker dB SL') + xlabel('gap') +% text(0.02,80,' TMC','backgroundColor','w') + end + title([num2str(BFs(BFno)) ' Hz']) + set(gca,'XTick',[ 0.1],'xTickLabel', { '0.1'}) +end + +% IFMCs +for BFno=1:length(compareBFs) + compareFreq=compareMaskerFreqs'*BFs(BFno); + subplot(2,1,2) + semilogx(compareFreq,compareIFMCs(:,BFno),'k:') + ylim([0 100]) + xlim([100 12000]) + grid on +end
--- a/multithreshold 1.46/printReport.m Fri Jul 01 12:59:47 2011 +0100 +++ b/multithreshold 1.46/printReport.m Fri Jul 08 13:48:27 2011 +0100 @@ -6,7 +6,7 @@ global experiment stimulusParameters betweenRuns withinRuns statsModel audio global LevittControl expGUIhandles - +global paramChanges global inputStimulusParams OMEParams DRNLParams global IHC_VResp_VivoParams IHCpreSynapseParams AN_IHCsynapseParams @@ -24,7 +24,9 @@ printReportGuide.showTracks=experiment.printTracks; printReportGuide.fileName=[]; if experiment.saveData - saveFileName=['savedData/' experiment.name '_' experiment.date '_' experiment.paradigm]; + saveFileName=... + ['savedData/' experiment.name '_' ... + experiment.date '_' experiment.paradigm]; else % save this data (just in case) saveFileName=['savedData/mostRecentResults']; @@ -84,6 +86,7 @@ msg=printTabTable(resultsTable, headers); addToMsg(msg,0) fprintf('\n') +disp(paramChanges) % sort tracks into the same order betweenRuns.levelTracks=betweenRuns.levelTracks(idx1); @@ -181,6 +184,7 @@ end fprintf('\nparadigm:\t%s\n ', experiment.paradigm) +disp(paramChanges) % ------------------------------------------------------- sortTablesForPrinting function table= sortTablesForPrinting(idx1,idx2, var1values,var2values, x)
--- a/multithreshold 1.46/subjGUI_MT.m Fri Jul 01 12:59:47 2011 +0100 +++ b/multithreshold 1.46/subjGUI_MT.m Fri Jul 08 13:48:27 2011 +0100 @@ -1349,11 +1349,11 @@ % without waiting for button press startNewRun(handles) - % show sample Rate on GUI; it must be set in MAPparams + % show sample Rate on GUI; it must be set in MAPparams ##?? set(expGUIhandles.textsampleRate,'string',... num2str(stimulusParameters.sampleRate)) - if experiment.singleShot + if experiment.singleShot % ##?? AN_IHCsynapseParams.plotSynapseContents=1; else AN_IHCsynapseParams.plotSynapseContents=0; @@ -1424,7 +1424,9 @@ function [modelResponse, MacGregorResponse]=MAPmodel global experiment stimulusParameters audio withinRuns -global outerMiddleEarParams DRNLParams AN_IHCsynapseParams +% global outerMiddleEarParams DRNLParams AN_IHCsynapseParams +global ICoutput ANdt dt savedBFlist ANprobRateOutput expGUIhandles +global paramChanges savePath=path; addpath(['..' filesep 'MAP'], ['..' filesep 'utilities']) @@ -1444,12 +1446,12 @@ % ---------------------------------------------- run Model MAPparamsName=experiment.name; showPlotsAndDetails=experiment.MAPplot; + +% important buried constant ##?? AN_spikesOrProbability='spikes'; AN_spikesOrProbability='probability'; % [response, method]=MAPsequenceSeg(audio, method, 1:8); -global ICoutput ANdt dt savedBFlist ANprobRateOutput expGUIhandles -global stimulusParameters experiment if sum(strcmp(experiment.ear,{'MAPmodelMultiCh', 'MAPmodelListen'})) % use BFlist specified in MAPparams file @@ -1458,7 +1460,8 @@ BFlist=stimulusParameters.targetFrequency; end paramChanges=get(expGUIhandles.editparamChanges,'string'); -eval(paramChanges) +% convert from string to a cell array +eval(paramChanges); MAP1_14(audio, stimulusParameters.sampleRate, BFlist,... MAPparamsName, AN_spikesOrProbability, paramChanges);
--- a/parameterStore/MAPparamsEndo.m Fri Jul 01 12:59:47 2011 +0100 +++ b/parameterStore/MAPparamsEndo.m Fri Jul 08 13:48:27 2011 +0100 @@ -129,7 +129,6 @@ % #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.Et= 0.05; % endocochlear potential (V) IHC_cilia_RPParams.Gk= 2e-008; % 1e-8 potassium conductance (S)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/parameterStore/MAPparamsJE.m Fri Jul 08 13:48:27 2011 +0100 @@ -0,0 +1,302 @@ +function method=MAPparamsJE ... + (BFlist, sampleRate, showParams, paramChanges) +% 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 IHC_cilia_RPParams +global IHC_VResp_VivoParams IHCpreSynapseParams AN_IHCsynapseParams +global MacGregorParams MacGregorMultiParams filteredSACFParams +global experiment % used by calls from multiThreshold only + + +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.006; % * N(all ICspikes) +% OMEParams.rateToAttenuationFactor=0; % * N(all ICspikes) + +% 'probability model': Ar based on AN firing probabilities (LSR) +OMEParams.rateToAttenuationFactorProb=0.01;% * N(all ANrates) +% OMEParams.rateToAttenuationFactorProb=0;% * 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=5e4; % DRNL.a=0 means no OHCs (no nonlinear path) +DRNLParams.a=2e4; % DRNL.a=0 means no OHCs (no nonlinear path) + +DRNLParams.b=8e-6; % *compression threshold raised compression +% DRNLParams.b=1; % b=1 means no compression + +DRNLParams.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=1000; % 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 = .01; % strength of MOC +% DRNLParams.rateToAttenuationFactor = 0; % strength of MOC +% 'probability' model: MOC based on AN spiking activity (HSR) +DRNLParams.rateToAttenuationFactorProb = .0055; % strength of MOC +% DRNLParams.rateToAttenuationFactorProb = .0; % strength of MOC +DRNLParams.MOCrateThresholdProb =70; % spikes/s probability only + +DRNLParams.MOCtau =.1; % smoothing for MOC + + +%% #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.03; % 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= 6e-9; % 2.5e-9 maximum conductance (Siemens) +IHC_cilia_RPParams.Ga= 1e-9; % 4.3e-9 fixed apical membrane conductance +IHC_cilia_RPParams.Ga= .8e-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.05; % 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=40e-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 + + +%% now accept last minute parameter changes required by the calling program +% paramChanges +if nargin>3 && ~isempty(paramChanges) + nChanges=length(paramChanges); + for idx=1:nChanges + eval(paramChanges{idx}) + end +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 + + % highlight parameter changes made locally + if nargin>3 && ~isempty(paramChanges) + fprintf('\n Local parameter changes:\n') + for i=1:length(paramChanges) + disp(paramChanges{i}) + end + end +end + +% for backward compatibility +experiment.comparisonData=[];
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/testAN.m Fri Jul 08 13:48:27 2011 +0100 @@ -0,0 +1,333 @@ +function vectorStrength=testAN(targetFrequency,BFlist, levels, ... + paramsName,paramChanges) +% testIHC used either for IHC I/O function ... +% or receptive field (doReceptiveFields=1) + +global IHC_VResp_VivoParams IHC_cilia_RPParams IHCpreSynapseParams +global AN_IHCsynapseParams + + global ANoutput ANdt CNoutput ICoutput ICmembraneOutput tauCas + global ARattenuation MOCattenuation + +dbstop if error +restorePath=path; + +addpath (['..' filesep 'MAP'], ['..' filesep 'utilities'], ... + ['..' filesep 'parameterStore'], ['..' filesep 'wavFileStore'],... + ['..' filesep 'testPrograms']) + +if nargin<5 + paramChanges=[]; +end + +if nargin<4 + paramsName='Normal'; +end + +if nargin<3 + levels=-10:10:80; + % levels=80:10:90; +end + +nLevels=length(levels); + +toneDuration=.2; +rampDuration=0.002; +silenceDuration=.02; +localPSTHbinwidth=0.001; + +% Use only the first frequency in the GUI targetFrequency box to defineBF +% targetFrequency=stimulusParameters.targetFrequency(1); +% BFlist=targetFrequency; + +AN_HSRonset=zeros(nLevels,1); +AN_HSRsaturated=zeros(nLevels,1); +AN_LSRonset=zeros(nLevels,1); +AN_LSRsaturated=zeros(nLevels,1); +CNLSRrate=zeros(nLevels,1); +CNHSRsaturated=zeros(nLevels,1); +ICHSRsaturated=zeros(nLevels,1); +ICLSRsaturated=zeros(nLevels,1); +vectorStrength=zeros(nLevels,1); + +AR=zeros(nLevels,1); +MOC=zeros(nLevels,1); + +% ANoutput=zeros(200,200); + +figure(15), clf +set(gcf,'position',[980 356 401 321]) +figure(5), clf +set(gcf,'position', [980 34 400 295]) +drawnow + +%% guarantee that the sample rate is at least 10 times the frequency +sampleRate=50000; +while sampleRate< 10* targetFrequency + sampleRate=sampleRate+10000; +end + +%% adjust sample rate so that the pure tone stimulus has an integer +%% numver of epochs in a period +dt=1/sampleRate; +period=1/targetFrequency; +ANspeedUpFactor=5; %anticipating MAP (needs clearing up) +ANdt=dt*ANspeedUpFactor; +ANdt=period/round(period/ANdt); +dt=ANdt/ANspeedUpFactor; +sampleRate=1/dt; +epochsPerPeriod=sampleRate*period; + +%% main computational loop (vary level) +levelNo=0; +for leveldB=levels + levelNo=levelNo+1; + + fprintf('%4.0f\t', leveldB) + amp=28e-6*10^(leveldB/20); + + time=dt:dt:toneDuration; + rampTime=dt:dt:rampDuration; + ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ... + ones(1,length(time)-length(rampTime))]; + ramp=ramp.*fliplr(ramp); + + silence=zeros(1,round(silenceDuration/dt)); + + % create signal (leveldB/ targetFrequency) + inputSignal=amp*sin(2*pi*targetFrequency'*time); + inputSignal= ramp.*inputSignal; + inputSignal=[silence inputSignal]; + + %% run the model + AN_spikesOrProbability='spikes'; + showPlotsAndDetails=0; + + + MAP1_14(inputSignal, 1/dt, BFlist, ... + paramsName, AN_spikesOrProbability, paramChanges); + + nTaus=length(tauCas); + + %LSR (same as HSR if no LSR fibers present) + [nANFibers nTimePoints]=size(ANoutput); + + numLSRfibers=nANFibers/nTaus; + numHSRfibers=numLSRfibers; + + LSRspikes=ANoutput(1:numLSRfibers,:); + PSTH=UTIL_PSTHmaker(LSRspikes, ANdt, localPSTHbinwidth); + PSTHLSR=mean(PSTH,1)/localPSTHbinwidth; % across fibers rates + PSTHtime=localPSTHbinwidth:localPSTHbinwidth:... + localPSTHbinwidth*length(PSTH); + AN_LSRonset(levelNo)= max(PSTHLSR); % peak in 5 ms window + AN_LSRsaturated(levelNo)= mean(PSTHLSR(round(length(PSTH)/2):end)); + + % HSR + HSRspikes= ANoutput(end- numHSRfibers+1:end, :); + PSTH=UTIL_PSTHmaker(HSRspikes, ANdt, localPSTHbinwidth); + PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only) + AN_HSRonset(levelNo)= max(PSTH); + AN_HSRsaturated(levelNo)= mean(PSTH(round(length(PSTH)/2): end)); + + figure(5), subplot(2,2,2) + hold off, bar(PSTHtime,PSTH, 'b') + hold on, bar(PSTHtime,PSTHLSR,'r') + ylim([0 1000]) + xlim([0 length(PSTH)*localPSTHbinwidth]) + set(gcf,'name',[num2str(BFlist), ' Hz: ' num2str(leveldB) ' dB']); + + % AN - CV + % CV is computed 5 times. Use the middle one (3) as most typical + cvANHSR= UTIL_CV(HSRspikes, ANdt); + + % AN - vector strength + PSTH=sum(HSRspikes); + [PH, binTimes]=UTIL_periodHistogram... + (PSTH, ANdt, targetFrequency); + VS=UTIL_vectorStrength(PH); + vectorStrength(levelNo)=VS; + disp(['sat rate= ' num2str(AN_HSRsaturated(levelNo)) ... + '; phase-locking VS = ' num2str(VS)]) + title(['AN HSR: CV=' num2str(cvANHSR(3),'%5.2f') ... + 'VS=' num2str(VS,'%5.2f')]) + + % CN - first-order neurons + + % CN LSR + [nCNneurons c]=size(CNoutput); + nLSRneurons=round(nCNneurons/nTaus); + CNLSRspikes=CNoutput(1:nLSRneurons,:); + PSTH=UTIL_PSTHmaker(CNLSRspikes, ANdt, localPSTHbinwidth); + PSTH=sum(PSTH)/nLSRneurons; + CNLSRrate(levelNo)=mean(PSTH(round(length(PSTH)/2):end))/localPSTHbinwidth; + + %CN HSR + MacGregorMultiHSRspikes=... + CNoutput(end-nLSRneurons:end,:); + PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, ANdt, localPSTHbinwidth); + PSTH=sum(PSTH)/nLSRneurons; + PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only) + + CNHSRsaturated(levelNo)=mean(PSTH(length(PSTH)/2:end)); + + figure(5), subplot(2,2,3) + bar(PSTHtime,PSTH) + ylim([0 1000]) + xlim([0 length(PSTH)*localPSTHbinwidth]) + cvMMHSR= UTIL_CV(MacGregorMultiHSRspikes, ANdt); + title(['CN CV= ' num2str(cvMMHSR(3),'%5.2f')]) + + % IC LSR + [nICneurons c]=size(ICoutput); + nLSRneurons=round(nICneurons/nTaus); + ICLSRspikes=ICoutput(1:nLSRneurons,:); + PSTH=UTIL_PSTHmaker(ICLSRspikes, ANdt, localPSTHbinwidth); + ICLSRsaturated(levelNo)=mean(PSTH(round(length(PSTH)/2):end))/localPSTHbinwidth; + + %IC HSR + MacGregorMultiHSRspikes=... + ICoutput(end-nLSRneurons:end,:); + PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, ANdt, localPSTHbinwidth); + PSTH=sum(PSTH)/nLSRneurons; + PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only) + + ICHSRsaturated(levelNo)=mean(PSTH(length(PSTH)/2:end)); + + AR(levelNo)=min(ARattenuation); + MOC(levelNo)=min(MOCattenuation(length(MOCattenuation)/2:end)); + + time=dt:dt:dt*size(ICmembraneOutput,2); + figure(5), subplot(2,2,4) + plot(time,ICmembraneOutput(2, 1:end),'k') + ylim([-0.07 0]) + xlim([0 max(time)]) + title(['IC ' num2str(leveldB,'%4.0f') 'dB']) + drawnow + + figure(5), subplot(2,2,1) + plot(20*log10(MOC), 'k'), + title(' MOC'), ylabel('dB attenuation') + ylim([-30 0]) + + +end % level +figure(5), subplot(2,2,1) +plot(levels,20*log10(MOC), 'k'), +title(' MOC'), ylabel('dB attenuation') +ylim([-30 0]) +xlim([0 max(levels)]) + +fprintf('\n') +toneDuration=2; +rampDuration=0.004; +silenceDuration=.02; +nRepeats=200; % no. of AN fibers +fprintf('toneDuration %6.3f\n', toneDuration) +fprintf(' %6.0f AN fibers (repeats)\n', nRepeats) +fprintf('levels') +fprintf('%6.2f\t', levels) +fprintf('\n') + + +% ---------------------------------------------------- display parameters + +figure(15), clf +nRows=2; nCols=2; + +% AN rate - level ONSET functions +subplot(nRows,nCols,1) +plot(levels,AN_LSRonset,'ro'), hold on +plot(levels,AN_HSRonset,'ko'), hold off +ylim([0 1000]), xlim([min(levels) max(levels)]) +ttl=['tauCa= ' num2str(IHCpreSynapseParams.tauCa)]; +title( ttl) +xlabel('level dB SPL'), ylabel('peak rate (sp/s)'), grid on +text(0, 800, 'AN onset', 'fontsize', 16) + +% AN rate - level ADAPTED function +subplot(nRows,nCols,2) +plot(levels,AN_LSRsaturated, 'ro'), hold on +plot(levels,AN_HSRsaturated, 'ko'), hold off +ylim([0 400]) +set(gca,'ytick',0:50:300) +xlim([min(levels) max(levels)]) +set(gca,'xtick',[levels(1):20:levels(end)]) +% grid on +ttl=[ 'spont=' num2str(mean(AN_HSRsaturated(1,:)),'%4.0f')... + ' sat=' num2str(mean(AN_HSRsaturated(end,1)),'%4.0f')]; +title( ttl) +xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)') +text(0, 340, 'AN adapted', 'fontsize', 16), grid on + +% CN rate - level ADAPTED function +subplot(nRows,nCols,3) +plot(levels,CNLSRrate, 'ro'), hold on +plot(levels,CNHSRsaturated, 'ko'), hold off +ylim([0 400]) +set(gca,'ytick',0:50:300) +xlim([min(levels) max(levels)]) +set(gca,'xtick',[levels(1):20:levels(end)]) +% grid on +ttl=[ 'spont=' num2str(mean(CNHSRsaturated(1,:)),'%4.0f') ' sat=' ... + num2str(mean(CNHSRsaturated(end,1)),'%4.0f')]; +title( ttl) +xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)') +text(0, 350, 'CN', 'fontsize', 16), grid on + +% IC rate - level ADAPTED function +subplot(nRows,nCols,4) +plot(levels,ICLSRsaturated, 'ro'), hold on +plot(levels,ICHSRsaturated, 'ko'), hold off +ylim([0 400]) +set(gca,'ytick',0:50:300) +xlim([min(levels) max(levels)]) +set(gca,'xtick',[levels(1):20:levels(end)]), grid on + + +ttl=['spont=' num2str(mean(ICHSRsaturated(1,:)),'%4.0f') ... + ' sat=' num2str(mean(ICHSRsaturated(end,1)),'%4.0f')]; +title( ttl) +xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)') +text(0, 350, 'IC', 'fontsize', 16) +set(gcf,'name',' AN CN IC rate/level') + +peakVectorStrength=max(vectorStrength); + +fprintf('\n') +disp('levels vectorStrength') +fprintf('%3.0f \t %6.4f \n', [levels; vectorStrength']) +fprintf('\n') +fprintf('Phase locking, max vector strength=\t %6.4f\n\n',... + max(vectorStrength)) + +allData=[ levels' AN_HSRonset AN_HSRsaturated... + AN_LSRonset AN_LSRsaturated ... + CNHSRsaturated CNLSRrate... + ICHSRsaturated ICLSRsaturated]; +fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR \tICLSR \n'); +UTIL_printTabTable(round(allData)) +fprintf('VS (phase locking)= \t%6.4f\n\n',... + max(vectorStrength)) + +UTIL_showStruct(IHC_cilia_RPParams, 'IHC_cilia_RPParams') +UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams') +UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams') + +fprintf('\n') +disp('levels vectorStrength') +fprintf('%3.0f \t %6.4f \n', [levels; vectorStrength']) +fprintf('\n') +fprintf('Phase locking, max vector strength= \t%6.4f\n\n',... + max(vectorStrength)) + +allData=[ levels' AN_HSRonset AN_HSRsaturated... + AN_LSRonset AN_LSRsaturated ... + CNHSRsaturated CNLSRrate... + ICHSRsaturated ICLSRsaturated]; +fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR \tICLSR \n'); +UTIL_printTabTable(round(allData)) +fprintf('VS (phase locking)= \t%6.4f\n\n',... + max(vectorStrength)) + +path(restorePath) \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/testANprob.m Fri Jul 08 13:48:27 2011 +0100 @@ -0,0 +1,199 @@ +function vectorStrength=testANprob(targetFrequency,BFlist, levels, ... + paramsName, paramChanges) +% testIHC used either for IHC I/O function ... +% or receptive field (doReceptiveFields=1) + +global IHC_VResp_VivoParams IHC_cilia_RPParams IHCpreSynapseParams +global AN_IHCsynapseParams + +global ANprobRateOutput dt tauCas +global ARattenuation MOCattenuation + +AN_spikesOrProbability='probability'; + +dbstop if error +addpath (['..' filesep 'MAP'], ['..' filesep 'utilities'], ... + ['..' filesep 'parameterStore'], ['..' filesep 'wavFileStore'],... + ['..' filesep 'testPrograms']) + +if nargin<5 + paramChanges=[]; +end + +if nargin<4 + paramsName='Normal'; +end + +if nargin<3 + levels=-10:10:80; +end + +nLevels=length(levels); + +toneDuration=.2; +rampDuration=0.002; +silenceDuration=.02; +localPSTHbinwidth=0.001; + +% Use only the first frequency in the GUI targetFrequency box to defineBF +% targetFrequency=stimulusParameters.targetFrequency(1); +% BFlist=targetFrequency; + +AN_HSRonset=zeros(nLevels,1); +AN_HSRsaturated=zeros(nLevels,1); +AN_LSRonset=zeros(nLevels,1); +AN_LSRsaturated=zeros(nLevels,1); + +AR=zeros(nLevels,1); +MOC=zeros(nLevels,1); + +figure(15), clf +set(gcf,'position',[607 17 368 321]) +drawnow + +%% guarantee that the sample rate is at least 10 times the frequency +sampleRate=50000; +while sampleRate< 10* targetFrequency + sampleRate=sampleRate+10000; +end + +%% adjust sample rate so that the pure tone stimulus has an integer +%% numver of epochs in a period +dt=1/sampleRate; +period=1/targetFrequency; + +%% main computational loop (vary level) +levelNo=0; +for leveldB=levels + levelNo=levelNo+1; + + fprintf('%4.0f\t', leveldB) + amp=28e-6*10^(leveldB/20); + + time=dt:dt:toneDuration; + rampTime=dt:dt:rampDuration; + ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ... + ones(1,length(time)-length(rampTime))]; + ramp=ramp.*fliplr(ramp); + + silence=zeros(1,round(silenceDuration/dt)); + + % create signal (leveldB/ targetFrequency) + inputSignal=amp*sin(2*pi*targetFrequency'*time); + inputSignal= ramp.*inputSignal; + inputSignal=[silence inputSignal]; + + %% run the model + showPlotsAndDetails=0; + + + MAP1_14(inputSignal, 1/dt, BFlist, ... + paramsName, AN_spikesOrProbability, paramChanges); + + nTaus=length(tauCas); + + %LSR (same as HSR if no LSR fibers present) + [nANFibers nTimePoints]=size(ANprobRateOutput); + + numLSRfibers=1; + numHSRfibers=numLSRfibers; + + LSRspikes=ANprobRateOutput(1:numLSRfibers,:); + PSTH=UTIL_PSTHmaker(LSRspikes, dt, localPSTHbinwidth); + PSTHLSR=PSTH/(localPSTHbinwidth/dt); % across fibers rates + PSTHtime=localPSTHbinwidth:localPSTHbinwidth:... + localPSTHbinwidth*length(PSTH); + AN_LSRonset(levelNo)= max(PSTHLSR); % peak in 5 ms window + AN_LSRsaturated(levelNo)= mean(PSTHLSR(round(length(PSTH)/2):end)); + + % HSR + HSRspikes= ANprobRateOutput(end- numHSRfibers+1:end, :); + PSTH=UTIL_PSTHmaker(HSRspikes, dt, localPSTHbinwidth); + PSTH=PSTH/(localPSTHbinwidth/dt); % sum across fibers (HSR only) + AN_HSRonset(levelNo)= max(PSTH); + AN_HSRsaturated(levelNo)= mean(PSTH(round(length(PSTH)/2): end)); + + figure(15), subplot(2,2,4) + hold off, bar(PSTHtime,PSTH, 'b') + hold on, bar(PSTHtime,PSTHLSR,'r') + ylim([0 1000]) + xlim([0 length(PSTH)*localPSTHbinwidth]) + set(gcf,'name',[num2str(BFlist), ' Hz: ' num2str(leveldB) ' dB']); + + AR(levelNo)=min(ARattenuation); + MOC(levelNo)=min(MOCattenuation(length(MOCattenuation)/2:end)); + + + figure(15), subplot(2,2,3) + plot(20*log10(MOC), 'k'), + title(' MOC'), ylabel('dB attenuation') + ylim([-30 0]) + + +end % level +figure(15), subplot(2,2,3) +plot(levels,20*log10(MOC), 'k'), +title(' MOC'), ylabel('dB attenuation') +ylim([-30 0]) +xlim([0 max(levels)]) + +fprintf('\n') +toneDuration=2; +rampDuration=0.004; +silenceDuration=.02; +nRepeats=200; % no. of AN fibers +fprintf('toneDuration %6.3f\n', toneDuration) +fprintf(' %6.0f AN fibers (repeats)\n', nRepeats) +fprintf('levels') +fprintf('%6.2f\t', levels) +fprintf('\n') + + +% ---------------------------------------------------- display parameters + + +nRows=2; nCols=2; + +% AN rate - level ONSET functions +subplot(nRows,nCols,1) +plot(levels,AN_LSRonset,'ro'), hold on +plot(levels,AN_HSRonset,'ko'), hold off +ylim([0 1000]), xlim([min(levels) max(levels)]) +ttl=['tauCa= ' num2str(IHCpreSynapseParams.tauCa)]; +title( ttl) +xlabel('level dB SPL'), ylabel('peak rate (sp/s)'), grid on +text(0, 800, 'AN onset', 'fontsize', 16) + +% AN rate - level ADAPTED function +subplot(nRows,nCols,2) +plot(levels,AN_LSRsaturated, 'ro'), hold on +plot(levels,AN_HSRsaturated, 'ko'), hold off +ylim([0 400]) +set(gca,'ytick',0:50:300) +xlim([min(levels) max(levels)]) +set(gca,'xtick',[levels(1):20:levels(end)]) +% grid on +ttl=[ 'spont=' num2str(mean(AN_HSRsaturated(1,:)),'%4.0f')... + ' sat=' num2str(mean(AN_HSRsaturated(end,1)),'%4.0f')]; +title( ttl) +xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)') +text(0, 340, 'AN adapted', 'fontsize', 16), grid on + +allData=[ levels' AN_HSRonset AN_HSRsaturated... + AN_LSRonset AN_LSRsaturated ]; +fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR \tICLSR \n'); +UTIL_printTabTable(round(allData)) + + +UTIL_showStruct(IHC_cilia_RPParams, 'IHC_cilia_RPParams') +UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams') +UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams') + +fprintf('\n') +disp('levels vectorStrength') + +allData=[ levels' AN_HSRonset AN_HSRsaturated... + AN_LSRonset AN_LSRsaturated ]; +fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR \tICLSR \n'); +UTIL_printTabTable(round(allData)) +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/testBM.m Fri Jul 08 13:48:27 2011 +0100 @@ -0,0 +1,164 @@ +function testBM (BMlocations, paramsName,... + relativeFrequencies, AN_spikesOrProbability, paramChanges) +% testBM generates input output functions for DRNL model for any number +% of locations. +% Computations are bast on a single channel model (channelBFs=BF) +% peak displacement (peakAmp) is measured. +% if DRNLParams.useMOC is chosen, the full model is run (slow) +% otherwise only DRNL is computed. +% Tuning curves are generated based on a range of frequencies reletove to +% the BF of the location. +% + +global DRNLParams + +if nargin<5 + paramChanges=[]; +end + +if nargin<4 + AN_spikesOrProbability='probability'; +end + +savePath=path; +addpath (['..' filesep 'utilities'],['..' filesep 'MAP']) + +levels=-10:10:90; nLevels=length(levels); +% levels= 50; nLevels=length(levels); + +% refBMdisplacement is the displacement of the BM at threshold +% 1 nm disp at threshold (9 kHz, Ruggero) +% ? adjust for frequency +refBMdisplacement= 1e-8; % adjusted for 10 nm at 1 kHz + +toneDuration=.200; +rampDuration=0.01; +silenceDuration=0.01; + +sampleRate=30000; + +dbstop if error +figure(3), clf +set(gcf,'position',[280 350 327 326]) +set(gcf,'name','DRNL - BM') +pause(0.1) + +finalSummary=[]; +nBFs=length(BMlocations); +BFno=0; plotCount=0; +for BF=BMlocations + BFno=BFno+1; + plotCount=plotCount+nBFs; + stimulusFrequencies=BF* relativeFrequencies; + nFrequencies=length(stimulusFrequencies); + + peakAmpBM=zeros(nLevels,nFrequencies); + peakAmpBMdB=NaN(nLevels,nFrequencies); + peakEfferent=NaN(nLevels,nFrequencies); + peakAREfferent=NaN(nLevels,nFrequencies); + + + levelNo=0; + for leveldB=levels + disp(['level= ' num2str(leveldB)]) + levelNo=levelNo+1; + + freqNo=0; + for frequency=stimulusFrequencies + freqNo=freqNo+1; + + % Generate stimuli + globalStimParams.FS=sampleRate; + globalStimParams.overallDuration=... + toneDuration+silenceDuration; % s + stim.phases='sin'; + stim.type='tone'; + stim.toneDuration=toneDuration; + stim.frequencies=frequency; + stim.amplitudesdB=leveldB; + stim.beginSilence=silenceDuration; + stim.rampOnDur=rampDuration; + % no offset ramp + stim.rampOffDur=rampDuration; + doPlot=0; + inputSignal=stimulusCreate(globalStimParams, stim, doPlot); + inputSignal=inputSignal(:,1)'; + + %% run the model + MAPparamsName=paramsName; + + global DRNLoutput MOCattenuation ARattenuation + MAP1_14(inputSignal, sampleRate, BF, ... + MAPparamsName, AN_spikesOrProbability, paramChanges); + + DRNLresponse=DRNLoutput; + peakAmp=max(max(... + DRNLresponse(:, end-round(length(DRNLresponse)/2):end))); + peakAmpBM(levelNo,freqNo)=peakAmp; + peakAmpBMdB(levelNo,freqNo)=... + 20*log10(peakAmp/refBMdisplacement); + peakEfferent(levelNo,freqNo)=min(min(MOCattenuation)); + peakAREfferent(levelNo,freqNo)=min(min(ARattenuation)); + + end % tone frequency + end % level + + %% analyses results and plot + + % BM I/O plot (top panel) + figure(3) + subplot(3,nBFs,BFno), cla + plot(levels,peakAmpBMdB, 'linewidth',2) + hold on, plot(levels, repmat(refBMdisplacement,1,length(levels))) + hold off + title(['BF=' num2str(BF,'%5.0f') ' - ' paramsName]) + xlabel('level') + % set(gca,'xtick',levels), grid on + if length(levels)>1,xlim([min(levels) max(levels)]), end + ylabel(['dB re:' num2str(refBMdisplacement,'%6.1e') 'm']) + ylim([-20 50]) + set(gca,'ytick',[-10 0 10 20 40]) + grid on + % legend({num2str(stimulusFrequencies')}, 'location', 'EastOutside') + UTIL_printTabTable([levels' peakAmpBMdB], ... + num2str([0 stimulusFrequencies]','%5.0f'), '%5.0f') + finalSummary=[finalSummary peakAmpBMdB]; + + % Tuning curve + if length(relativeFrequencies)>2 + figure(3), subplot(3,nBFs, 2*nBFs+BFno) + % contour(stimulusFrequencies,levels,peakAmpBM,... + % [refBMdisplacement refBMdisplacement],'r') + contour(stimulusFrequencies,levels,peakAmpBM,... + refBMdisplacement.*[1 5 10 50 100]) + title(['tuning curve at ' num2str(refBMdisplacement) 'm']); + ylabel('level (dB) at reference') + xlim([100 10000]) + grid on + hold on + set(gca,'xscale','log') + end + + + % MOC contribution + figure(3) + subplot(3,nBFs,nBFs+BFno), cla + plot(levels,20*log10(peakEfferent), 'linewidth',2) + ylabel('MOC (dB attenuation)'), xlabel('level') + title(['peak MOC: model= ' AN_spikesOrProbability]) + grid on + if length(levels)>1, xlim([min(levels) max(levels)]), end + + % AR contribution + hold on + plot(levels,20*log10(peakAREfferent), 'r') + hold off + +end % best frequency + +UTIL_showStructureSummary(DRNLParams, 'DRNLParams', 10) + + UTIL_printTabTable([levels' finalSummary], ... + num2str([0 stimulusFrequencies]','%5.0f'), '%5.0f') + +path(savePath);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/testFM.m Fri Jul 08 13:48:27 2011 +0100 @@ -0,0 +1,255 @@ +function testFM(BFlist,paramsName,showPSTHs,paramChanges) +% specify whether you want AN 'probability' or 'spikes' +% spikes is more realistic but takes longer +% refractory effect is included only for spikes +% + +% specify the AN ANresponse bin width (normally 1 ms). +% This can influence the measure of the AN onset rate based on the +% largest bin in the ANresponse +% +% Demonstration is based on Harris and Dallos + +global inputStimulusParams outerMiddleEarParams DRNLParams +global IHC_VResp_VivoParams IHCpreSynapseParams AN_IHCsynapseParams + +dbstop if error + +restorePath=path; +addpath (['..' filesep 'MAP'], ['..' filesep 'utilities'], ... + ['..' filesep 'parameterStore'], ['..' filesep 'wavFileStore'],... + ['..' filesep 'testPrograms']) + +if nargin<3 + paramChanges=[]; +end + +% masker and probe levels are relative to this threshold +thresholdAtCF=10; % dB SPL +thresholdAtCF=5; % dB SPL + +if nargin<1, showPSTHs=1;end + +sampleRate=50000; + +% fetch BF from GUI: use only the first target frequency +maskerFrequency=BFlist; +maskerDuration=.1; + +targetFrequency=maskerFrequency; +probeLeveldB=20+thresholdAtCF; % H&D use 20 dB SL/ TMC uses 10 dB SL +probeDuration=0.008; % HD use 15 ms probe (fig 3). +probeDuration=0.015; % HD use 15 ms probe (fig 3). + +rampDuration=.001; % HD use 1 ms linear ramp +initialSilenceDuration=0.02; +finalSilenceDuration=0.05; % useful for showing recovery + +maskerLevels=[-80 10 20 30 40 60 ] + thresholdAtCF; +% maskerLevels=[-80 40 60 ] + thresholdAtCF; + +dt=1/sampleRate; + +figure(7), clf +set(gcf,'position',[613 36 360 247]) +set(gcf,'name', ['forward masking: thresholdAtCF=' num2str(thresholdAtCF)]) + +if showPSTHs + figure(8), clf + set(gcf,'name', 'Harris and Dallos simulation') + set(gcf,'position',[980 36 380 249]) +end + +% Plot Harris and Dallos result for comparison +gapDurations=[0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.3]; +HDmaskerLevels=[+10 +20 +30 +40 +60]; +HDresponse=[ + 1 1 1 1 1 1 1 1; + 0.8 0.82 0.81 0.83 0.87 0.95 1 1; + 0.48 0.5 0.54 0.58 0.7 0.85 0.95 1; + 0.3 0.31 0.35 0.4 0.5 0.68 0.82 0.94; + 0.2 0.27 0.27 0.29 0.42 0.64 0.75 0.92; + 0.15 0.17 0.18 0.23 0.3 0.5 0.6 0.82]; + +figure(7) +semilogx(gapDurations,HDresponse,'o'), hold on +legend(strvcat(num2str(maskerLevels')),'location','southeast') +title([ 'masker dB: ' num2str(HDmaskerLevels)]) + +%% Run the trials +maxProbeResponse=0; +levelNo=0; +resultsMatrix=zeros(length(maskerLevels), length(gapDurations)); +summaryFiringRates=[]; + +% initial silence +time=dt: dt: initialSilenceDuration; +initialSilence=zeros(1,length(time)); + +% probe +time=dt: dt: probeDuration; +amp=28e-6*10^(probeLeveldB/20); +probe=amp*sin(2*pi.*targetFrequency*time); +% ramps +% catch rampTime error +if rampDuration>0.5*probeDuration, rampDuration=probeDuration/2; end +rampTime=dt:dt:rampDuration; +% raised cosine ramp +ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ... + ones(1,length(time)-length(rampTime))]; +% onset ramp +probe=probe.*ramp; +% offset ramp makes complete ramp for probe +ramp=fliplr(ramp); +% apply ramp mask to probe. Probe does not change below +probe=probe.*ramp; + +% final silence +time=dt: dt: finalSilenceDuration; +finalSilence=zeros(1,length(time)); + +PSTHplotCount=0; +longestSignalDuration=initialSilenceDuration + maskerDuration +... + max(gapDurations) + probeDuration + finalSilenceDuration ; +for maskerLeveldB=maskerLevels + levelNo=levelNo+1; + allDurations=[]; + allFiringRates=[]; + + %masker + time=dt: dt: maskerDuration; + masker=sin(2*pi.*maskerFrequency*time); + % masker ramps + if rampDuration>0.5*maskerDuration + % catch ramp duration error + rampDuration=maskerDuration/2; + end + % NB masker ramp (not probe ramp) + rampTime=dt:dt:rampDuration; + % raised cosine ramp + ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi))... + ones(1,length(time)-length(rampTime))]; + % onset ramp + masker=masker.*ramp; + % offset ramp + ramp=fliplr(ramp); + % apply ramp + masker=masker.*ramp; + amp=28e-6*10^(maskerLeveldB/20); + maskerPa=amp*masker; + + for gapDuration=gapDurations + time=dt: dt: gapDuration; + gap=zeros(1,length(time)); + + inputSignal=... + [initialSilence maskerPa gap probe finalSilence]; + + % ********************************** run MAP model + + global ANprobRateOutput tauCas ... + + showPlotsAndDetails=0; + +AN_spikesOrProbability='probability'; + +MAP1_14(inputSignal, 1/dt, targetFrequency, ... + paramsName, AN_spikesOrProbability, paramChanges); + + [nFibers c]=size(ANprobRateOutput); + nLSRfibers=nFibers/length(tauCas); + ANresponse=ANprobRateOutput(end-nLSRfibers:end,:); + ANresponse=sum(ANresponse)/nLSRfibers; + + % analyse results + probeStart=initialSilenceDuration+maskerDuration+gapDuration; + PSTHbinWidth=dt; + responseDelay=0.005; + probeTimes=probeStart+responseDelay:... + PSTHbinWidth:probeStart+probeDuration+responseDelay; + probeIDX=round(probeTimes/PSTHbinWidth); + probePSTH=ANresponse(probeIDX); + firingRate=mean(probePSTH); + % NB this only works if you start with the lowest level masker + maxProbeResponse=max([maxProbeResponse firingRate]); + allDurations=[allDurations gapDuration]; + allFiringRates=[allFiringRates firingRate]; + + %% show PSTHs + if showPSTHs + nLevels=length(maskerLevels); + nDurations=length(gapDurations); + figure(8) + PSTHbinWidth=1e-3; + PSTH=UTIL_PSTHmaker(ANresponse, dt, PSTHbinWidth); + PSTH=PSTH*dt/PSTHbinWidth; + PSTHplotCount=PSTHplotCount+1; + subplot(nLevels,nDurations,PSTHplotCount) + probeTime=PSTHbinWidth:PSTHbinWidth:... + PSTHbinWidth*length(PSTH); + bar(probeTime, PSTH) + if strcmp(AN_spikesOrProbability, 'spikes') + ylim([0 500]) + else + ylim([0 500]) + end + xlim([0 longestSignalDuration]) + grid on + + if PSTHplotCount< (nLevels-1) * nDurations+1 + set(gca,'xticklabel',[]) + end + + if ~isequal(mod(PSTHplotCount,nDurations),1) + set(gca,'yticklabel',[]) + else + ylabel([num2str(maskerLevels... + (round(PSTHplotCount/nDurations) +1))]) + end + + if PSTHplotCount<=nDurations + title([num2str(1000*gapDurations(PSTHplotCount)) 'ms']) + end + end % showPSTHs + + end % gapDurations duration + summaryFiringRates=[summaryFiringRates allFiringRates']; + + figure(7), hold on + semilogx(allDurations, allFiringRates/maxProbeResponse) + ylim([0 1]), hold on + ylim([0 inf]), xlim([min(gapDurations) max(gapDurations)]) + xlim([0.001 1]) + pause(0.1) % to allow for CTRL/C interrupts + resultsMatrix(levelNo,:)=allFiringRates; +end % maskerLevel + +disp('delay/ rates') +disp(num2str(round( [1000*allDurations' summaryFiringRates]))) + +% replot with best adjustment +% figure(34), clf% use for separate plot +figure(7), clf +peakProbe=max(max(resultsMatrix)); +resultsMatrix=resultsMatrix/peakProbe; +semilogx(gapDurations,HDresponse,'o'), hold on +title(['FrMsk: probe ' num2str(probeLeveldB)... + 'dB SL: peakProbe=' num2str(peakProbe,'%5.0f') ' sp/s']) +xlabel('gap duration (s)'), ylabel ('probe response') +semilogx(allDurations, resultsMatrix'), ylim([0 1]) +ylim([0 inf]), +xlim([0.001 1]) +legend(strvcat(num2str(maskerLevels'-thresholdAtCF)), -1) + +% ------------------------------------------------- display parameters +disp(['parameter file was: ' paramsName]) +fprintf('\n') +% UTIL_showStruct(inputStimulusParams, 'inputStimulusParams') +% UTIL_showStruct(outerMiddleEarParams, 'outerMiddleEarParams') +% UTIL_showStruct(DRNLParams, 'DRNLParams') +% UTIL_showStruct(IHC_VResp_VivoParams, 'IHC_VResp_VivoParams') +UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams') +UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams') + + +path(restorePath); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/testOME.m Fri Jul 08 13:48:27 2011 +0100 @@ -0,0 +1,90 @@ +function testOME(paramsName, paramChanges) + +savePath=path; +addpath (['..' filesep 'utilities'],['..' filesep 'MAP']) + +if nargin<2 + paramChanges=[]; +end + +sampleRate=50000; + +dt=1/sampleRate; +leveldBSPL=80; % dB SPL as used by Huber (may trigger AR) +amp=10^(leveldBSPL/20)*28e-6; +duration=.05; +time=dt: dt: duration; + +%% Comparison data (human) +% These data are taken directly from Huber 2001 (Fig. 4) +HuberFrequencies=[600 800 1000 2000 3000 4000 6000 8000]; +HuberDisplacementAt80dBSPL=[1.5E-9 1.5E-09 1.5E-09 1.0E-09 7.0E-10 ... + 3.0E-10 2.0E-10 1.0E-10]; % m; +% HuberVelocityAt80dBSPL= 2*pi*HuberFrequencies.*HuberDisplacementAt80dBSPL; + +figure(2), clf, subplot(2,1,1) +set(2,'position',[5 349 268 327]) +semilogx(HuberFrequencies, 20*log10(HuberDisplacementAt80dBSPL/1e-10),... + 'ko', 'MarkerFaceColor','k', 'Marker','o', 'markerSize',6) +hold on + +%% Generate test stimulus ................................................................. + +% independent test using discrete frequencies +peakResponses=[]; +peakTMpressure=[]; +frequencies=[200 400 HuberFrequencies 10000]; +for toneFrequency=frequencies + inputSignal=amp*sin(2*pi*toneFrequency*time); + + showPlotsAndDetails=0; + AN_spikesOrProbability='probability'; + % switch off AR & MOC (Huber's patients were deaf) + paramChanges{1}='OMEParams.rateToAttenuationFactorProb=0;'; + paramChanges{2}='DRNLParams.rateToAttenuationFactorProb = 0;'; + + global OMEoutput OMEextEarPressure TMoutput ARattenuation + % BF is irrelevant + MAP1_14(inputSignal, sampleRate, -1, ... + paramsName, AN_spikesOrProbability, paramChanges); + + peakDisplacement=max(OMEoutput(floor(end/2):end)); + peakResponses=[peakResponses peakDisplacement]; + + peakTMpressure=[peakTMpressure max(OMEextEarPressure)]; + disp([' AR attenuation (dB): ' num2str(20*log10(min(ARattenuation)))]) +end + +%% Report +disp('frequency displacement(m)') +% disp(num2str([frequencies' peakResponses'])) +fprintf('%6.0f \t%10.3e\n',[frequencies' peakResponses']') + +% 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 +title(['stapes at ' num2str(leveldBSPL) ' (NB deaf)']) +ylabel('disp: dB re 1e-10m') +xlabel('stimulus frequency (Hz)') +legend({'Huber et al','model'},'location','southWest') +set(gcf,'name','OME') + +% external ear resonance +figure(2), subplot(2,1,2),hold off +semilogx(frequencies, 20*log10(peakTMpressure/28e-6)-leveldBSPL,... + 'k', 'linewidth',2) +xlim([100 10000]) %, ylim([-10 30]) +grid on +title(['External ear resonances' ]) +ylabel('dB') +xlabel('frequency') +set(gcf,'name','OME: external resonances') +% ---------------------------------------------------------- display parameters +disp(['parameter file was: ' paramsName]) +fprintf('\n') + +path(savePath);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/testPeriphery.m Fri Jul 08 13:48:27 2011 +0100 @@ -0,0 +1,161 @@ +function testPeriphery (BMlocations, paramsName) +% testBM generates input output functions for DRNL model for any number +% of locations. +% Computations are bast on a single channel model (channelBFs=BF) +% peak displacement (peakAmp) is measured. +% if DRNLParams.useMOC is chosen, the full model is run (slow) +% otherwise only DRNL is computed. +% Tuning curves are generated based on a range of frequencies reletove to +% the BF of the location. +% + +global DRNLParams + +savePath=path; + +addpath (['..' filesep 'utilities'],['..' filesep 'MAP']) + +levels=-10:10:90; nLevels=length(levels); +% levels= 50; nLevels=length(levels); + +% relativeFrequencies=[0.25 .5 .75 1 1.25 1.5 2]; +relativeFrequencies=1; + +% refBMdisplacement is the displacement of the BM at threshold +% 1 nm disp at threshold (9 kHz, Ruggero) +refBMdisplacement= 1e-8; % adjusted for 10 nm at 1 kHz + +toneDuration=.200; +rampDuration=0.01; +silenceDuration=0.01; + +sampleRate=30000; + +dbstop if error +figure(3), clf +% set(gcf,'position',[276 33 331 645]) +set(gcf,'name','DRNL - BM') + +finalSummary=[]; +nBFs=length(BMlocations); +BFno=0; plotCount=0; +for BF=BMlocations + BFno=BFno+1; + plotCount=plotCount+nBFs; + stimulusFrequencies=BF* relativeFrequencies; + nFrequencies=length(stimulusFrequencies); + + peakAmpBM=zeros(nLevels,nFrequencies); + peakAmpBMdB=NaN(nLevels,nFrequencies); + peakEfferent=NaN(nLevels,nFrequencies); + peakAREfferent=NaN(nLevels,nFrequencies); + + + levelNo=0; + for leveldB=levels + disp(['level= ' num2str(leveldB)]) + levelNo=levelNo+1; + + freqNo=0; + for frequency=stimulusFrequencies + freqNo=freqNo+1; + + % Generate stimuli + globalStimParams.FS=sampleRate; + globalStimParams.overallDuration=... + toneDuration+silenceDuration; % s + stim.type='tone'; + stim.phases='sin'; + stim.toneDuration=toneDuration; + stim.frequencies=frequency; + stim.amplitudesdB=leveldB; + stim.beginSilence=silenceDuration; + stim.rampOnDur=rampDuration; + stim.rampOffDur=rampDuration; + doPlot=0; + inputSignal=stimulusCreate(globalStimParams, stim, doPlot); + inputSignal=inputSignal(:,1)'; + + %% run the model + MAPparamsName=paramsName; + AN_spikesOrProbability='probability'; + % spikes are slow but can be used to study MOC using IC units + % AN_spikesOrProbability='spikes'; + + global DRNLoutput MOCattenuation ARattenuation IHCoutput + MAP1_14(inputSignal, sampleRate, BF, ... + MAPparamsName, AN_spikesOrProbability); + + DRNLresponse=IHCoutput; + peakAmp=max(max(... + DRNLresponse(:, end-round(length(DRNLresponse)/2):end))); + peakAmpBM(levelNo,freqNo)=peakAmp; + if peakAmp>0 + peakAmpBMdB(levelNo,freqNo)=... + 20*log10(peakAmp/refBMdisplacement); + else + peakAmpBMdB(levelNo,freqNo)=peakAmp; + end + peakEfferent(levelNo,freqNo)=min(min(MOCattenuation)); + peakAREfferent(levelNo,freqNo)=min(min(ARattenuation)); + + end % tone frequency + end % level + + %% analyses results and plot + + % BM I/O plot (top panel) + figure(3) + subplot(3,nBFs,BFno), cla + plot(levels,peakAmpBMdB, 'linewidth',2) + hold on, plot(levels, repmat(refBMdisplacement,1,length(levels))) + hold off + title(['BF=' num2str(BF,'%5.0f') ' - ' paramsName]) + xlabel('level') + % set(gca,'xtick',levels), grid on + if length(levels)>1,xlim([min(levels) max(levels)]), end + ylabel(['dB re:' num2str(refBMdisplacement,'%6.1e') 'm']) + ylim([-20 50]) + set(gca,'ytick',[-10 0 10 20 40]) + % legend({num2str(stimulusFrequencies')}, 'location', 'EastOutside') + UTIL_printTabTable([levels' peakAmpBMdB], ... + num2str([0 stimulusFrequencies]','%5.0f'), '%5.0f') + finalSummary=[finalSummary peakAmpBMdB]; + + % Tuning curve + if length(relativeFrequencies)>2 + figure(3), subplot(3,nBFs, nBFs+BFno) + % contour(stimulusFrequencies,levels,peakAmpBM,... + % [refBMdisplacement refBMdisplacement],'r') + contour(stimulusFrequencies,levels,peakAmpBM,... + refBMdisplacement.*[1 5 10 50 100]) + title(['tuning curve at ' num2str(refBMdisplacement) 'm']); + ylabel('level (dB) at reference') + xlim([100 10000]) + hold on + set(gca,'xscale','log') + end + + + % MOC contribution + figure(3) + subplot(3,nBFs,2*nBFs+BFno), cla + plot(levels,20*log10(peakEfferent), 'linewidth',2) + ylabel('MOC (dB attenuation)'), xlabel('level') + title(['peak MOC: model= ' AN_spikesOrProbability]) + grid on + if length(levels)>1, xlim([min(levels) max(levels)]), end + + % AR contribution + hold on + plot(levels,20*log10(peakAREfferent), 'r') + hold off + +end % best frequency + +UTIL_showStructureSummary(DRNLParams, 'DRNLParams', 10) + + UTIL_printTabTable([levels' finalSummary], ... + num2str([0 stimulusFrequencies]','%5.0f'), '%5.0f') + +path(savePath); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/testPhaseLocking.m Fri Jul 08 13:48:27 2011 +0100 @@ -0,0 +1,51 @@ +function testPhaseLocking(paramsName, paramChanges) + +if nargin<2 + paramChanges=[]; +end + +testFrequencies=[250 500 1000 2000 4000 8000]; +levels=50:10:80; +figure(14), clf +set(gcf,'position', [980 36 383 321]) +set(gcf,'name', 'phase locking') +allStrengths=zeros(length(testFrequencies), length(levels)); +peakVectorStrength=zeros(1,length(testFrequencies)); +freqCount=0; +for targetFrequency=testFrequencies; + %single test + freqCount=freqCount+1; + vectorStrength=... + testAN(targetFrequency,targetFrequency, levels,... + paramsName, paramChanges); + allStrengths(freqCount,:)=vectorStrength'; + peakVectorStrength(freqCount)=max(vectorStrength'); +end +%% plot results +figure(14) +subplot(2,1,2) +plot(levels,allStrengths) +xlabel('levels') +ylabel('vector strength') +legend (num2str(testFrequencies'),'location','eastOutside') + +subplot(2,1,1) +semilogx(testFrequencies,peakVectorStrength) +grid on +title ('peak vector strength') +xlabel('frequency') +ylabel('vector strength') + +johnson=[250 0.79 +500 0.82 +1000 0.8 +2000 0.7 +4000 0.25 +5500 0.05 +]; +hold on +plot(johnson(:,1),johnson(:,2),'o') +legend({'model','Johnson 80'},'location','eastOutside') +hold off + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/testPhysiology.m Fri Jul 08 13:48:27 2011 +0100 @@ -0,0 +1,21 @@ +function testPhysiology(BF,paramsName, paramChanges) + +restorePath=path; +addpath (['..' filesep 'MAP']) + +if nargin<3 + paramChanges=[]; +end + +testOME(paramsName,paramChanges) +relativeFrequencies=[0.25 .5 .75 1 1.25 1.5 2]; +testBM (BF, paramsName,relativeFrequencies,'spikes', paramChanges) +testRP(BF,paramsName,paramChanges) +testSynapse(BF,paramsName,paramChanges) +testFM(BF,paramsName,1,paramChanges) +testPhaseLocking(paramsName,paramChanges) +testAN(BF,BF, -10:10:80,paramsName,paramChanges); +figure(14) +figure(15) + +path(restorePath)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/testPhysiologyProb.m Fri Jul 08 13:48:27 2011 +0100 @@ -0,0 +1,19 @@ +function testPhysiologyProb(BF,paramsName, paramChanges) + +restorePath=path; +addpath (['..' filesep 'MAP']) + +if nargin<3 + paramChanges=[]; +end + +testOME(paramsName, paramChanges) +relativeFrequencies=[0.25 .5 .75 1 1.25 1.5 2]; +testBM (BF, paramsName,relativeFrequencies,'probability', paramChanges) +testRP(BF,paramsName, paramChanges) +testSynapse(BF,paramsName, paramChanges) +testANprob(BF,BF, -10:10:80,paramsName, paramChanges); + +figure(4) + +path(restorePath)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/testRF.m Fri Jul 08 13:48:27 2011 +0100 @@ -0,0 +1,295 @@ +function testRF +% testIHC used either for IHC I/O function or receptive field (doReceptiveFields=1) + +global experiment method stimulusParameters expGUIhandles +global inputStimulusParams IHC_ciliaParams +global IHC_VResp_VivoParams IHCpreSynapseParams AN_IHCsynapseParams +dbstop if error +% set(expGUIhandles.pushbuttonStop, 'backgroundColor', [.941 .941 .941]) + +addpath (['..' filesep 'MAP'], ['..' filesep 'utilities'], ... + ['..' filesep 'parameterStore'], ['..' filesep 'wavFileStore'],... + ['..' filesep 'testPrograms']) + +targetFrequency=stimulusParameters.targetFrequency(1); + +sampleRate=50000; +doReceptiveFields=1; + +toneDuration=.05; +rampDuration=0.004; +silenceDuration=.02; + +nRepeats=100; % no. of AN fibers + +plotGraphsForIHC=1; +% number of MacGregor units is set in the parameter file. + +if doReceptiveFields + % show all receptive field + frequencies=targetFrequency* [ 0.5 0.7 0.9 1 1.1 1.3 1.6]; + levels=0:20:80; nLevels=length(levels); + figure(14), clf + figure(15), clf +else + % show only I/O function at BF + frequencies=targetFrequency; + levels=-20:10:90; + % levels=10:.25:13; + % levels=-20:1:-15 + nLevels=length(levels); +% figure(13), clf, +% set (gcf, 'name', ['IHC/AN input/output' num2str(AN_IHCsynapseParams.numFibers) ' repeats']) +% drawnow +end +nFrequencies=length(frequencies); + +IHC_RP_peak=zeros(nLevels,nFrequencies); +IHC_RP_min=zeros(nLevels,nFrequencies); +IHC_RP_dc=zeros(nLevels,nFrequencies); +AN_HSRonset=zeros(nLevels,nFrequencies); +AN_HSRsaturated=zeros(nLevels,nFrequencies); +AN_LSRonset=zeros(nLevels,nFrequencies); +AN_LSRsaturated=zeros(nLevels,nFrequencies); +CNLSRsaturated=zeros(nLevels,nFrequencies); +CNHSRsaturated=zeros(nLevels,nFrequencies); +ICHSRsaturated=zeros(nLevels,nFrequencies); +ICLSRsaturated=zeros(nLevels,nFrequencies); + + +levelNo=0; PSTHplotCount=0; +for leveldB=levels + fprintf('%4.0f\t', leveldB) + levelNo=levelNo+1; + amp=28e-6*10^(leveldB/20); + + freqNo=0; + for frequency=frequencies + + paramFunctionName=['method=MAPparams' experiment.name ... + '(' num2str(targetFrequency) ');' ]; + eval(paramFunctionName); % read parameters afresh each pass + + dt=method.dt; + time=dt:dt:toneDuration; + rampTime=dt:dt:rampDuration; + ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ones(1,length(time)-length(rampTime))]; + ramp=ramp.*fliplr(ramp); + + silence=zeros(1,round(silenceDuration/dt)); + + toneStartptr=length(silence)+1; + toneMidptr=toneStartptr+round(toneDuration/(2*dt)) -1; + toneEndptr=toneStartptr+round(toneDuration/dt) -1; + + % create signal (leveldB/ frequency) + freqNo=freqNo+1; + inputSignal=amp*sin(2*pi*frequency'*time); + inputSignal= ramp.*inputSignal; + inputSignal=[silence inputSignal silence]; + + if doReceptiveFields % receptive field + method.plotGraphs= 0; % plot only PSTHs + else + method.plotGraphs= plotGraphsForIHC; % show progress + end + + targetChannelNo=1; + + % force parameters + % the number of AN fibers at each BF + AN_IHCsynapseParams.numFibers= nRepeats; + AN_IHCsynapseParams. mode= 'spikes'; + AN_IHCsynapseParams.plotSynapseContents=0; + AN_IHCsynapseParams.PSTHbinWidth=.001; + + method.DRNLSave=1; + method.IHC_cilia_RPSave=1; + method.PSTHbinWidth=1e-3; % useful 1-ms default for all PSTHs + method.AN_IHCsynapseSave=1; + method.MacGregorMultiSave=1; + method.MacGregorSave=1; + method.dt=dt; + + moduleSequence=[1:8]; + + global ANdt ARAttenuation TMoutput OMEoutput ... + DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV... + IHCoutput ANprobRateOutput ANoutput savePavailable tauCas ... + CNoutput ICoutput ICmembraneOutput ICfiberTypeRates MOCattenuation + +AN_spikesOrProbability='spikes'; +AN_spikesOrProbability='probability'; +MAPparamsName='Normal'; + +MAP1_14(inputSignal, 1/dt, targetFrequency, ... + MAPparamsName, AN_spikesOrProbability); + + % RP + IHC_RPData=IHC_cilia_output; + IHC_RPData=IHCoutput(targetChannelNo,:); + IHC_RP_peak(levelNo,freqNo)=max(IHC_RPData(toneStartptr:toneEndptr)); + IHC_RP_min(levelNo,freqNo)=min(IHC_RPData(toneStartptr:toneEndptr)); + IHC_RP_dc(levelNo,freqNo)=mean(IHC_RPData(toneStartptr:toneEndptr)); + + % AN next + AN_IHCsynapseAllData=ANoutput; + method.PSTHbinWidth=0.001; + + nTaus=length(tauCas); + numANfibers=size(ANoutput,1); + numLSRfibers=numANfibers/nTaus; + + %LSR (same as HSR if no LSR fibers present) + channelPtr1=(targetChannelNo-1)*numANfibers+1; + channelPtr2=channelPtr1+numANfibers-1; + LSRspikes=AN_IHCsynapseAllData(channelPtr1:channelPtr2,:); + method.dt=method.AN_IHCsynapsedt; + PSTH=UTIL_PSTHmaker(LSRspikes, method); + PSTH=sum(PSTH,1); % sum across fibers (HSR only) + PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1; + PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1; + PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1; + AN_LSRonset(levelNo,freqNo)=max(max(PSTH))/(method.PSTHbinWidth*method.numANfibers); + AN_LSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr))/(method.numANfibers*toneDuration/2); + + % HSR + channelPtr1=numLSRfibers+(targetChannelNo-1)*method.numANfibers+1; + channelPtr2=channelPtr1+method.numANfibers-1; + HSRspikes=AN_IHCsynapseAllData(channelPtr1:channelPtr2,:); + method.dt=method.AN_IHCsynapsedt; + PSTH=UTIL_PSTHmaker(HSRspikes, method); + PSTH=sum(PSTH,1); % sum across fibers (HSR only) + PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1; + PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1; + PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1; + AN_HSRonset(levelNo,freqNo)=max(max(PSTH))/(method.PSTHbinWidth*method.numANfibers); + AN_HSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr))/(method.numANfibers*toneDuration/2); + [cvANHSR, cvTimes, allTimeStamps, allISIs]= UTIL_CV(HSRspikes, method.AN_IHCsynapsedt); + + PSTHplotCount=PSTHplotCount+1; + if doReceptiveFields % receptive field for HSR only + figure(14), set(gcf,'name','AN') + plotReceptiveFields(method, PSTH, PSTHplotCount, levels, frequencies) + ylim([0 method.numANfibers]) + xlabel(['CV= ' num2str(max(cvANHSR),'%4.2f')],'fontsize',8) + end % doReceptiveFields + + % CN + MacGregorMultiAllData=method.MacGregorMultiData; + numLSRfibers=method.McGMultinNeuronsPerBF*length(method.nonlinCF)* (nTaus-1); + + %LSR (same as HSR if no LSR fibers present) + channelPtr1=(targetChannelNo-1)*method.McGMultinNeuronsPerBF+1; + channelPtr2=channelPtr1+method.McGMultinNeuronsPerBF-1; + MacGregorMultiLSRspikes=MacGregorMultiAllData(channelPtr1:channelPtr2,:); + method.dt=method.MacGregorMultidt; + PSTH=UTIL_PSTHmaker(MacGregorMultiLSRspikes, method); + PSTH=sum(PSTH,1); % sum across fibers (HSR only) + PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1; + PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1; + PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1; + CNLSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr)); + CNLSRsaturated(levelNo,freqNo)=CNLSRsaturated(levelNo,freqNo)... + /((toneDuration/2)*method.McGMultinNeuronsPerBF); + + %HSR + channelPtr1=numLSRfibers+(targetChannelNo-1)*method.McGMultinNeuronsPerBF+1; + channelPtr2=channelPtr1+method.McGMultinNeuronsPerBF-1; + MacGregorMultiHSRspikes=MacGregorMultiAllData(channelPtr1:channelPtr2,:); + method.dt=method.MacGregorMultidt; + PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, method); + PSTH=sum(PSTH,1); % sum across fibers (HSR only) + PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1; + PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1; + PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1; + CNHSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr)); + CNHSRsaturated(levelNo,freqNo)=CNHSRsaturated(levelNo,freqNo)... + /((toneDuration/2)*method.McGMultinNeuronsPerBF); + [cvMMHSR, cvTimes, allTimeStamps, allISIs]= UTIL_CV(MacGregorMultiHSRspikes, method.MacGregorMultidt); + + if doReceptiveFields % receptive field + figure(15), set(gcf,'name','CN HSR input') + plotReceptiveFields(method, PSTH, PSTHplotCount, levels, frequencies) + ylim([0 method.McGMultinNeuronsPerBF]) + xlabel(['CV= ' num2str(max(cvMMHSR),'%4.2f')],'fontsize',8) + end + + MacGregorAllData=method.MacGregorData; + numLSRfibers=length(method.nonlinCF)* (nTaus-1); + + %LSR (same as HSR if no LSR fibers present) + channelPtr1=targetChannelNo; + MacGregorLSR=MacGregorAllData(channelPtr1,:); + method.dt=method.MacGregordt; + PSTH=UTIL_PSTHmaker(MacGregorLSR, method); + PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1; + PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1; + PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1; + ICLSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr)); + ICLSRsaturated(levelNo,freqNo)=ICLSRsaturated(levelNo,freqNo)/(toneDuration/2); + + %LSR (same as HSR if no LSR fibers present) + channelPtr1=numLSRfibers+targetChannelNo; + MacGregorHSR=MacGregorAllData(channelPtr1,:); + method.dt=method.MacGregordt; + PSTH=UTIL_PSTHmaker(MacGregorHSR, method); + PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1; + PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1; + PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1; + ICHSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr)); + ICHSRsaturated(levelNo,freqNo)=ICHSRsaturated(levelNo,freqNo)/(toneDuration/2); + [cvICHSR, cvTimes, allTimeStamps, allISIs]= UTIL_CV(MacGregorHSR, method.MacGregordt); + +% if doReceptiveFields % receptive field +% figure(16), set(gcf,'name','IC HSR input') +% plotReceptiveFields(method, PSTH, PSTHplotCount, levels, frequencies) +% ylim([0 method.McGMultinNeuronsPerBF]) +% xlabel(['CV= ' num2str(max(cvICHSR),'%4.2f')],'fontsize',8) +% end + end % frequency +end % level +fprintf('\n') +toneDuration=2; +rampDuration=0.004; +silenceDuration=.02; +nRepeats=200; % no. of AN fibers +fprintf('toneDuration %6.3f\n', toneDuration) +fprintf(' %6.0f AN fibers (repeats)\n', nRepeats) +fprintf('levels') +fprintf('%6.2f\t', levels) +fprintf('\n') + + +% ---------------------------------------------------------- display parameters +disp(['parameter file was: ' experiment.name]) +fprintf('\n') +UTIL_showStruct(IHC_VResp_VivoParams, 'IHC_cilia_RPParams') +UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams') +UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams') + + + +function plotReceptiveFields(method, PSTH, PSTHplotCount, levels, frequencies) + +% show PSTH for each level/frequency combination +nLevels=length(levels); +nFrequencies=length(frequencies); + +PSTHtime=method.PSTHbinWidth:method.PSTHbinWidth:method.PSTHbinWidth*length(PSTH); +subplot(nLevels,nFrequencies,PSTHplotCount) +bar(PSTHtime, PSTH) +xlim([0 max(PSTHtime)]) +% write axis labels only at left and bottom +if PSTHplotCount< (nLevels-1) * nFrequencies+1 + set(gca,'xticklabel',[]) +end +if ~isequal(mod(PSTHplotCount,nFrequencies),1) + set(gca,'yticklabel',[]) +else + ylabel([num2str(levels(round(PSTHplotCount/nFrequencies) +1)) ' dB']) +end +% add titles only on top row +if PSTHplotCount<=nFrequencies + title([num2str(frequencies(PSTHplotCount)) ' Hz']) +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/testRP.m Fri Jul 08 13:48:27 2011 +0100 @@ -0,0 +1,205 @@ +function testRP(BFs,MAPparamsName,paramChanges) +% testIHC used for IHC I/O function + +global experiment method inputStimulusParams +global stimulusParameters IHC_VResp_VivoParams IHC_cilia_RPParams +savePath=path; +addpath (['..' filesep 'utilities'],['..' filesep 'MAP']) +dbstop if error + +figure(4), clf, +set (gcf, 'name', ['IHC']) +set(gcf,'position',[613 354 360 322]) +drawColors='rgbkmcy'; +drawnow + +if nargin<3 + paramChanges=[]; +end + +levels=-20:10:100; +nLevels=length(levels); +toneDuration=.05; +silenceDuration=.01; +sampleRate=50000; +dt=1/sampleRate; + +allIHC_RP_peak=[]; +allIHC_RP_dc=[]; + +for BFno=1:length(BFs) + BF=BFs(BFno); + targetFrequency=BF; + % OR + %Patuzzi and Sellick test (see ELP & AEM, 2006) + % targetFrequency=100; + + IHC_RP_peak=zeros(nLevels,1); + IHC_RP_min=zeros(nLevels,1); + IHC_RP_dc=zeros(nLevels,1); + + time=dt:dt:toneDuration; + + rampDuration=0.004; + rampTime=dt:dt:rampDuration; + ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ... + ones(1,length(time)-length(rampTime))]; + ramp=ramp.*fliplr(ramp); + + silence=zeros(1,round(silenceDuration/dt)); + + toneStartptr=length(silence)+1; + toneMidptr=toneStartptr+round(toneDuration/(2*dt)) -1; + toneEndptr=toneStartptr+round(toneDuration/dt) -1; + + levelNo=0; + for leveldB=levels + levelNo=levelNo+1; + % replicate at all levels + amp=28e-6*10^(leveldB/20); + + %% create signal (leveldB/ frequency) + inputSignal=amp*sin(2*pi*targetFrequency'*time); + inputSignal= ramp.*inputSignal; + inputSignal=[silence inputSignal silence]; + inputStimulusParams.sampleRate=1/dt; +% global IHC_ciliaParams + + %% disable efferent for fast processing + method.DRNLSave=1; + method.IHC_cilia_RPSave=1; + method.IHCpreSynapseSave=1; + method.IHC_cilia_RPSave=1; + method.segmentDuration=-1; + moduleSequence=1:4; + + %% run the model + global DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV... + IHCoutput + AN_spikesOrProbability='probability'; + + MAP1_14(inputSignal, sampleRate, BF, ... + MAPparamsName, AN_spikesOrProbability, paramChanges); + + % DRNL + DRNLoutput=DRNLoutput; + DRNL_peak(levelNo,1)=max(DRNLoutput(toneMidptr:toneEndptr)); + DRNL_min(levelNo,1)=min(DRNLoutput(toneMidptr:toneEndptr)); + DRNL_dc(levelNo,1)=mean(DRNLoutput(toneMidptr:toneEndptr)); + + % cilia + IHC_ciliaData=IHC_cilia_output; + IHC_ciliaData=IHC_ciliaData; + IHC_cilia_peak(levelNo,1)=... + max(IHC_ciliaData(toneMidptr:toneEndptr)); + IHC_cilia_min(levelNo,1)=... + min(IHC_ciliaData(toneMidptr:toneEndptr)); + IHC_cilia_dc(levelNo,1)=... + mean(IHC_ciliaData(toneMidptr:toneEndptr)); + + % RP + IHC_RPData=IHCoutput; + IHC_RPData=IHC_RPData; + IHC_RP_peak(levelNo,1)=... + max(IHC_RPData(toneMidptr:toneEndptr)); + IHC_RP_min(levelNo,1)=... + min(IHC_RPData(toneMidptr:toneEndptr)); + IHC_RP_dc(levelNo,1)=... + mean(IHC_RPData(toneMidptr:toneEndptr)); + + end % level + + + disp(['parameter file was: ' MAPparamsName]) + fprintf('\n') + + %% plot DRNL + subplot(2,2,1) +% referenceDisp= 9e-9*1000/targetFrequency; +% plot(levels,20*log10(DRNL_peak/referenceDisp), drawColors(BFno), ... +% 'linewidth',2), hold on + referenceDisp=10e-9; + plot(levels,20*log10(DRNL_peak/referenceDisp), drawColors(BFno), ... + 'linewidth',2), hold on + title([' DRNL peak: ' num2str(BFs) ' Hz']) + ylabel ('log10DRNL(m)'), xlabel('dB SPL') + xlim([min(levels) max(levels)]), ylim([-10 50]) + grid on + + %% plot cilia displacement + figure(4) + subplot(2,2,2) + restingIHC_cilia=IHCrestingCiliaCond; + plot(levels, IHC_cilia_peak,'k', 'linewidth',2), hold on + plot(levels, IHC_cilia_min,'r', 'linewidth',2) + hold on, + plot([min(levels) max(levels)], ... + [restingIHC_cilia restingIHC_cilia], 'g') + title(' IHC apical cond.') + ylabel ('IHCcilia(conductance)'), xlabel('dB SPL') + xlim([min(levels) max(levels)]) + grid on + + %% plot receptor potentials + figure(4) + subplot(2,2,3) + % RP I/O function min and max + restingRP=IHC_RP_peak(1); + toPlot= [fliplr(IHC_RP_min(:,1)') IHC_RP_peak(:,1)']; + microPa= 28e-6*10.^(levels/20); + microPa=[-fliplr(microPa) microPa]; + plot(microPa,toPlot, drawColors(BFno), 'linewidth',2) + % ylim([0 300]) + + %% Dallos and Harris data + dallosx=[-0.9 -0.1 -0.001 0.001 0.01 0.9]; + dallosy=[-8 -7.8 -6.5 11 16.5 22]/1000 + restingRP; + hold on, plot(dallosx,dallosy, 'o') + plot([-1 1], [restingRP restingRP], 'r') + title(' Dallos(86) data at 800 Hz') + ylabel ('receptor potential(V)'), xlabel('Pa') + ylim([-0.08 -0.02]), xlim([-1 1]) + grid on + + %% RP I/O function min and max + figure(4) + subplot(2,2,4) + restingRP=IHC_RP_peak(1); + peakRP=max(IHC_RP_peak); + plot(levels, IHC_RP_peak,drawColors(BFno), 'linewidth',2) + hold on + plot(levels, IHC_RP_dc, [drawColors(BFno) ':'], 'linewidth',2) + hold on, + plot([min(levels) max(levels)], [restingRP restingRP], 'r') + xlim([min(levels) max(levels)]) + % animal data + sndLevel=[5 15 25 35 45 55 65 75]; +RPanimal=restingRP+[0.5 2 4.6 5.8 6.4 7.2 8 10.2]/1000; +% could be misleading when restingRP changes +RPanimal=-0.060+[0.5 2 4.6 5.8 6.4 7.2 8 10.2]/1000; +hold on, plot(sndLevel,RPanimal,'o') + + grid on + title(['Et= ' num2str(IHC_cilia_RPParams.Et) ': RP data 7 kHz Patuzzi']) + ylabel ('RP(V)'), xlabel('dB SPL') + ylim([-0.08 -0.04]) + allIHC_RP_peak=[allIHC_RP_peak IHC_RP_peak]; + allIHC_RP_dc=[allIHC_RP_dc IHC_RP_dc]; + + fprintf('level\t peak\t DC\n') + UTIL_printTabTable([levels' IHC_RP_peak IHC_RP_dc]) + % disp(['restingIHC_cilia= ' num2str(restingIHC_cilia)]) + fprintf('peakRP= \t%6.3f', peakRP) + fprintf('\nrestingRP= \t%6.3f', restingRP) + fprintf('\ndifference= \t%6.3f\n', (peakRP-restingRP)) + drawnow +end +% UTIL_showStruct(IHC_VResp_VivoParams, 'IHC_VResp_VivoParams') +UTIL_showStruct(IHC_cilia_RPParams, 'IHC_cilia_RPParams') + fprintf('level\t peak\n') + UTIL_printTabTable([levels' allIHC_RP_peak]) + fprintf('level\t DC\n') + UTIL_printTabTable([levels' allIHC_RP_dc]) + +path(savePath); +disp(paramChanges)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/testSynapse.m Fri Jul 08 13:48:27 2011 +0100 @@ -0,0 +1,87 @@ +function testSynapse(BFlist,paramsName, paramChanges) +% testSynapse tracks the quantity of available transmitter vesicles +% the computations are single channel using the first frequency +% in the targetFrequency box of the expGUI. +% For, speed this function uses only probability and HSR fibers. +% This cannot be changed because of the way AN_IHCsynapse is coded.This). + +global experiment method IHCpreSynapseParams +global AN_IHCsynapseParams stimulusParameters +savePath=path; +addpath (['..' filesep 'utilities'],['..' filesep 'MAP']) + +if nargin<3 + paramChanges=[]; +end + +figure(6),clf +plotColors='rbgkrbgkrbgkrbgkrbgkrbgk'; +set(gcf,'position',[5 32 264 243]) + +sampleRate=5e4; dt=1/sampleRate; + +maskerLevels=-0:10:100; + +targetFrequency=BFlist; + +silenceDuration=0.015; +maskerDuration=0.1; +recoveryDuration=0.15; +rampDuration=0.004; + +maskerTime=dt:dt:maskerDuration; + +rampTime=dt:dt:rampDuration; +ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ... + ones(1,length(maskerTime)-length(rampTime))]; +ramp=ramp.*fliplr(ramp); + +initialSilence=zeros(1,round(silenceDuration/dt)); +recoverySilence=zeros(1,round(recoveryDuration/dt)); + +signal=sin(2*pi*targetFrequency'*maskerTime); +signal= ramp.*signal; +signal=[initialSilence signal recoverySilence]; + +levelCount=0; +qtMatrix=[]; +for leveldB=maskerLevels + levelCount=levelCount+1; + + amp=28e-6*10^(leveldB/20); + inputSignal=amp*signal; + + AN_spikesOrProbability='probability'; + showPlotsAndDetails=0; + + global savePavailable + + MAP1_14(inputSignal, 1/dt, targetFrequency, ... + paramsName, AN_spikesOrProbability, paramChanges); + + % ignore LSR channels (if any) at the top of the matrix + qt=savePavailable(end, :); + + synapsedt=dt; + time=synapsedt:synapsedt:synapsedt*length(qt); + + figure(6) + qtMatrix=[qtMatrix; qt]; + plot(time,qt, plotColors(levelCount)) + hold on + xlim([0 max(time)]) + ylim([0 AN_IHCsynapseParams.M]) +end + +set(gcf,'name','pre-synaptic available transmitter') +title(['q - available vesicles:' num2str(BFlist) ' Hz']) +legend(strvcat(num2str(maskerLevels')),'location','southeast') +legend boxoff +grid on + +figure(88), [c,H]=contour(time, maskerLevels,qtMatrix); clabel(c, H); +set(gcf,'position',[ 276 31 328 246]) +xlabel('time'), ylabel('maskerLevels') +title('contour plot of available transmitter') + +path(savePath);