comparison userProgramsRM/test_Dolan_and_Nuttall.m @ 38:c2204b18f4a2 tip

End nov big change

+function test_Dolan_and_Nuttall
+% test_MAP1_14 is a general purpose test routine that can be adjusted to
+% test a number of different applications of MAP1_14
+%
+% A range of options are supplied in the early part of the program
+%
+% One use of the function is to create demonstrations; filenames <demoxx>
+%  to illustrate particular features
+%
+% #1
+% Identify the file (in 'MAPparamsName') containing the model parameters
+%
+% #2
+% Identify the kind of model required (in 'AN_spikesOrProbability').
+%  A full brainstem model (spikes) can be computed or a shorter model
+%  (probability) that computes only so far as the auditory nerve
+%
+% #3
+% Choose between a tone signal or file input (in 'signalType')
+%
+% #4
+% Set the signal rms level (in leveldBSPL)
+%
+% #5
+% Identify the channels in terms of their best frequencies in the vector
+%  BFlist.
+%
+% Last minute changes to the parameters fetched earlier can be made using
+%  the cell array of strings 'paramChanges'.
+%  Each string must have the same format as the corresponding line in the
+%  file identified in 'MAPparamsName'
+%
+% When the demonstration is satisfactory, freeze it by renaming it <demoxx>
+
+global dt dtSpikes  savedBFlist saveAN_spikesOrProbability saveMAPparamsName...
+    savedInputSignal OMEextEarPressure TMoutput OMEoutput ARattenuation ...
+    DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
+    IHCoutput ANprobRateOutput ANoutput savePavailable ANtauCas  ...
+    CNtauGk CNoutput  ICoutput ICmembraneOutput ICfiberTypeRates ...
+    MOCattenuation
+global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams
+global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams
+global ICrate
+
+
+dbstop if error
+restorePath=path;
+addpath (['..' filesep 'MAP'],    ['..' filesep 'wavFileStore'], ...
+    ['..' filesep 'utilities'])
+
+%%  #1 parameter file name
+MAPparamsName='Normal';
+
+
+%% #2 probability (fast) or spikes (slow) representation
+AN_spikesOrProbability='spikes';
+%   or
+% AN_spikesOrProbability='probability';
+
+
+%% #3 pure tone, harmonic sequence or speech file input
+signalType= 'tones';
+toneFrequency= 4000;            % or a pure tone (Hz)
+
+sampleRate= 44100;          % must agree with noise
+duration=0.010;                 % seconds
+beginSilence=0.010;
+endSilence=0.010;
+rampDuration=.001;              % raised cosine ramp (seconds)
+noiseRampDuration=0.002;
+
+%   or
+% harmonic sequence (Hz)
+% F0=210;
+% toneFrequency= F0:F0:8000;
+
+%   or
+% signalType= 'file';
+% fileName='twister_44kHz';
+
+
+
+% %% #4 rms level
+% % signal details
+% leveldBSPL= 80;                  % dB SPL (80 for Lieberman)
+% leveldBSPLNoise=-30;
+
+%% #5 number of channels in the model
+%   21-channel model (log spacing)
+numChannels=21;
+lowestBF=250; 	highestBF= 8000;
+BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels));
+
+% %   or specify your own channel BFs
+% numChannels=1;
+% BFlist=toneFrequency;
+
+
+%% #6 change model parameters
+
+paramChanges={};
+
+% Parameter changes can be used to change one or more model parameters
+%  *after* the MAPparams file has been read
+% This example declares only one fiber type with a calcium clearance time
+% constant of 80e-6 s (HSR fiber) when the probability option is selected.
+% paramChanges={'AN_IHCsynapseParams.ANspeedUpFactor=5;', ...
+%     'IHCpreSynapseParams.tauCa=86e-6; '};
+% paramChanges={'DRNLParams.MOCtauProb =.25;', ...
+%     'DRNLParams.rateToAttenuationFactorProb = 0.02; '};
+
+paramChanges={'AN_IHCsynapseParams.numFibers=	50; ',...
+    'DRNLParams.MOCtauProb =.15;', ...
+    'DRNLParams.rateToAttenuationFactorProb = 0.00; '};
+
+% paramChanges={'AN_IHCsynapseParams.numFibers=	50; ',...
+% 'DRNLParams.rateToAttenuationFactorProb = -0.007;'};
+
+
+%% delare 'showMap' options to control graphical output
+showMapOptions.printModelParameters=1;   % prints all parameters
+showMapOptions.showModelOutput=0;       % plot of all stages
+showMapOptions.printFiringRates=1;      % prints stage activity levels
+showMapOptions.showACF=0;               % shows SACF (probability only)
+showMapOptions.showEfferent=1;          % tracks of AR and MOC
+showMapOptions.surfProbability=1;       % 2D plot of HSR response
+showMapOptions.surfSpikes=1;            % 2D plot of spikes histogram
+showMapOptions.ICrates=0;               % IC rates by CNtauGk
+
+% disable certain silly options
+if strcmp(AN_spikesOrProbability, 'spikes')
+    % avoid nonsensical options
+    showMapOptions.surfProbability=0;
+    showMapOptions.showACF=0;
+end
+
+if strcmp(signalType, 'file')
+    % needed for labeling plot
+    showMapOptions.fileName=fileName;
+else
+    showMapOptions.fileName=[];
+end
+
+fprintf('\n')
+disp([num2str(numChannels) ' channel model: ' AN_spikesOrProbability])
+disp('Computing ...')
+
+%%systematic
+probeLevels=30:10:80;
+noiseLevels=[-100 30];
+noRepeats=10;
+
+% probeLevels=80;
+% noiseLevels=[-30];
+% noRepeats=10;
+
+peakCAPs=zeros(4,length(probeLevels));
+
+for noiseCondition=1:length(noiseLevels)
+    leveldBSPLNoise=noiseLevels(noiseCondition);
+    levelNo=0;
+    for probeLevel=probeLevels
+        leveldBSPL=probeLevel;
+        levelNo=levelNo+1;
+        summedCAP=[];
+        for repeatNo= 1:noRepeats
+            disp(['repeat no: ' num2str(repeatNo)])
+            %% Generate stimuli
+
+            switch signalType
+                case 'tones'
+                    % Create pure tone stimulus
+                    dt=1/sampleRate; % seconds
+                    time=dt: dt: duration;
+                    inputSignal=sum(sin(2*pi*toneFrequency'*time), 1);
+                    amp=10^(leveldBSPL/20)*28e-6;   % converts to Pascals (peak)
+                    inputSignal=amp*inputSignal;
+                    % apply ramps
+                    % catch rampTime error
+                    if rampDuration>0.5*duration, rampDuration=duration/2; end
+                    rampTime=dt:dt:rampDuration;
+                    ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
+                        ones(1,length(time)-length(rampTime))];
+                    inputSignal=inputSignal.*ramp;
+                    ramp=fliplr(ramp);
+                    inputSignal=inputSignal.*ramp;
+                    % add silence
+                    intialSilence= zeros(1,round(beginSilence/dt));
+                    finalSilence= zeros(1,round(endSilence/dt));
+                    inputSignal= [intialSilence inputSignal finalSilence];
+
+                    %         [inputNoise sampleRateN]=wavread('babble');
+                    [inputNoise sampleRateN]=wavread('white noise');
+                    inputNoise=inputNoise(1:length(inputSignal));
+                    inputNoise=inputNoise(:,1);
+                    targetRMS=20e-6*10^(leveldBSPLNoise/20);
+                    rms=(mean(inputNoise.^2))^0.5;
+                    amp=targetRMS/rms;
+                    inputNoise=inputNoise*amp;
+                    time=dt: dt: dt*length(inputNoise);
+                    rampTime=dt:dt:noiseRampDuration;
+                    ramp=[0.5*(1+cos(2*pi*rampTime/(2*noiseRampDuration)+pi)) ...
+                        ones(1,length(time)-length(rampTime))];
+                    inputNoise=inputNoise'.*ramp;
+                    ramp=fliplr(ramp);
+                    inputNoise=inputNoise.*ramp;
+
+                    inputSignal=inputSignal+inputNoise;
+                    intialSilence= zeros(1,round(beginSilence/dt));
+                    finalSilence= zeros(1,round(endSilence/dt));
+                    inputSignal= [intialSilence inputSignal finalSilence];
+
+                    toneOnset=2*beginSilence;
+
+                    figure(2), subplot(3,1,1)
+                    time=dt:dt:dt*length(inputSignal);
+                    plot(time,inputSignal,'k')
+
+                case 'file'
+                    %% file input simple or mixed
+                    [inputSignal sampleRate]=wavread(fileName);
+                    dt=1/sampleRate;
+                    inputSignal=inputSignal(:,1);
+                    targetRMS=20e-6*10^(leveldBSPL/20);
+                    rms=(mean(inputSignal.^2))^0.5;
+                    amp=targetRMS/rms;
+                    inputSignal=inputSignal*amp;
+                    intialSilence= zeros(1,round(0.1/dt));
+                    finalSilence= zeros(1,round(0.2/dt));
+                    inputSignal= [intialSilence inputSignal' finalSilence];
+
+            end
+
+
+            %% run the model
+            tic
+
+            MAP1_14(inputSignal, sampleRate, BFlist, ...
+                MAPparamsName, AN_spikesOrProbability, paramChanges);
+
+
+            %% the model run is now complete. Now display the results
+                %                 UTIL_showMAP(showMapOptions, paramChanges)
+
+                wholeNerveCAP  = UTIL_CAPgenerator...
+                    (ANoutput, dtSpikes, BFlist, AN_IHCsynapseParams.numFibers, 1);
+
+                if isempty(summedCAP)
+                    summedCAP=wholeNerveCAP;
+                else
+                    summedCAP=summedCAP+wholeNerveCAP;
+                end
+
+                switch AN_spikesOrProbability
+                    case 'spikes'
+                        ANoutput = sum(ANoutput, 1);
+                    case 'probability'
+                        ANoutput = ANprobRateOutput(13+21,:);
+                end
+                figure(2), subplot(3,1,2), plot(ANoutput)
+                spikeTimes=dtSpikes:dtSpikes:dtSpikes* length(wholeNerveCAP);
+                figure(2), subplot(3,1,3), plot(spikeTimes,summedCAP/repeatNo)
+                ylim([-50 50])
+        end % repeat
+
+        spikeTimes=dtSpikes:dtSpikes:dtSpikes* length(wholeNerveCAP);
+        idx=find(spikeTimes>toneOnset & ...
+            spikeTimes>toneOnset+duration+.005);
+        averageCAP=summedCAP/repeatNo;
+        peakCAP=max(averageCAP(idx));
+        peakCAPs(noiseCondition,levelNo)=peakCAPs(noiseCondition,levelNo)+ peakCAP;
+
+        if strcmp(signalType,'tones')
+            disp(['duration=' num2str(duration)])
+            disp(['level=' num2str(leveldBSPL)])
+            disp(['toneFrequency=' num2str(toneFrequency)])
+            disp(['leveldBSPLNoise=' num2str(leveldBSPLNoise)])
+
+            disp(['attenuation factor =' ...
+                num2str(DRNLParams.rateToAttenuationFactor, '%5.3f') ])
+            disp(['attenuation factor (probability)=' ...
+                num2str(DRNLParams.rateToAttenuationFactorProb, '%5.3f') ])
+            disp(AN_spikesOrProbability)
+        end
+
+
+        disp([ 'peak CAP ' num2str(peakCAP)])
+
+        for i=1:length(paramChanges)
+            disp(paramChanges{i})
+        end
+    end % probe level
+    figure(9), subplot(3,1,3), plot(probeLevels,peakCAPs)
+end % condition
+%%
+
+path(restorePath)
+