comparison testPrograms/testLibermanMOC_DPOAE.m @ 38:c2204b18f4a2 tip

End nov big change

+function testLibermanMOC_DPOAE
+
+% compares MOC response to LIberman's 1996 data for DPOAE reduction with
+%  contralateral tone stimulation.
+% This program is used mainly as a check on the time constants involved.
+%
+% NB very different time constants are required for 'spikes' and
+% 'probability'
+
+% global dt ANdt  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 dt dtSpikes   saveAN_spikesOrProbability ANprobRateOutput ICoutput    
+
+global DRNLParams
+
+LibermanData=[
+    2	0.2;
+2.1	0.19;2.2	0.18;2.3	0.18;2.4	0.16;2.5	0.15;2.6	0.15;2.7 0.15;
+2.8	0.12;2.9	0.12;3	0.1;3.1	0.1;3.2	0.05;3.3	0.05;3.4	0;3.5	-0.1;
+3.6	-0.4;3.7	-1.2;3.8	-1.6;3.9	-1.8;4	-1.85;4.1	-2;4.2	-2.05;
+4.3	-2.05;4.4	-2.15;4.5	-2.2;4.6	-2.15;4.7	-2.1;4.8	-2.15;4.9 -2.2;
+5	-2.1;5.1	-2.1;5.2	-2.25;5.3	-2.1;5.4	-2.15;5.5	-2.1;5.6 -2.15;
+5.7	-2.1;5.8	-2.2;5.9	-2.05;6	-2.15;6.1	-2.05;6.2 -2;6.3 -2.2;6.4 -2.1;
+6.5	-2.05;6.6	-2.05;6.7	-2.05;6.8 -2.2;6.9 -2.1;7	-2.05;7.1 -2.05;7.2	-0.7;
+7.3	-0.1;7.4	0;7.5	0.1;7.6	0.2;7.7	0.35;7.8	0.2;7.9	0.15;8	0.2;8.1	0.15;8.2	0.15;
+8.3	0.15;8.4	0.12;8.5	0.1;8.6	0.09;8.7	0.08;8.8	0.07;8.9	0.06;9	0.05;
+];
+
+% Backus2006Data: time	bilateral	contralateral	ipsilateral
+% all % max microPascals
+Backus2006Data=[
+    100	20	15	10;
+    200	36	24	14;
+    300	44	30	18;
+    400	46	32	20;
+    500	48	34	22;
+    1000	50	36	24;
+    1500	52	37	25;
+    2000	54	38	27
+    ];
+
+steadyMinimum=mean(LibermanData(LibermanData(:,1)>4 & LibermanData(:,1)<7,2));
+
+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';
+sampleRate= 50000;
+rampDuration=.005;              % raised cosine ramp (seconds)
+toneFrequency= 1000;            % or a pure tone (Hz)
+duration=3.6;                   % Lieberman test
+beginSilence=1;                 % 1 for Lieberman
+endSilence=1;                   % 1 for Lieberman
+
+%% #4 rms level
+% signal details
+leveldBSPL= 80;                  % dB SPL (80 for Lieberman)
+
+
+%% #5 number of channels in the model
+
+numChannels=1;
+BFlist=toneFrequency;
+
+
+%% #6 change model parameters
+paramChanges={};
+
+%% delare 'showMap' options to control graphical output
+showMapOptions.printModelParameters=1;   % prints all parameters
+showMapOptions.showModelOutput=1;       % 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.surfAN=1;                % 2D plot of HSR response 
+showMapOptions.surfSpikes=0;            % 2D plot of spikes histogram
+showMapOptions.ICrates=0;               % IC rates by CNtauGk
+
+%% Generate stimuli
+
+        % 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];
+
+%% run the model
+tic
+fprintf('\n')
+disp(['Signal duration= ' num2str(length(inputSignal)/sampleRate)])
+disp([num2str(numChannels) ' channel model: ' AN_spikesOrProbability])
+disp('Computing ...')
+
+MAP1_14(inputSignal, sampleRate, BFlist, ...
+    MAPparamsName, AN_spikesOrProbability, paramChanges);
+
+
+%% the model run is now complete. Now display the results
+UTIL_showMAP(showMapOptions)
+
+if strcmp(signalType,'tones')
+    disp(['duration=' num2str(duration)])
+    disp(['level=' num2str(leveldBSPL)])
+    disp(['toneFrequency=' num2str(toneFrequency)])
+    disp(['attenuation factor =' ...
+        num2str(DRNLParams.rateToAttenuationFactor, '%5.3f') ])
+    disp(['attenuation factor (probability)=' ...
+        num2str(DRNLParams.rateToAttenuationFactorProb, '%5.3f') ])
+    disp(AN_spikesOrProbability)
+end
+disp(paramChanges)
+
+
+
+%% superimpose Lieberman (1996) data
+
+global  MOCattenuation
+MOCdB=20*log10(MOCattenuation);
+MOCtime=dt:dt:dt*length(MOCdB);
+
+% scale up DPOAE results to the running maximum MOC dB
+steadyMOCminimum=mean(MOCdB(MOCtime>2 & MOCtime<4.5));
+scalar=steadyMOCminimum/steadyMinimum;
+
+figure(90), clf
+plot(MOCtime,MOCdB), hold on
+plot(LibermanData(:,1)-2.5,scalar*LibermanData(:,2),'r:','linewidth',4)
+legend({'MAP', 'DPOAE'},'location', 'east')
+title('Compare Liberman (1996) DPOAE data with MAP MOC')
+xlabel('time (s)'), ylabel('MOC attenuation/ DPOAE reduction')
+if strcmp(saveAN_spikesOrProbability,'probability')
+text(0,2,['MOCtau= ' num2str(DRNLParams.MOCtauProb)])
+else
+text(0,2,['MOCtau= ' num2str(DRNLParams.MOCtau)])
+end
+set(gcf, 'name', 'Liberman compare')
+
+PSTHbinwidth=0.001;
+
+% show the source of the MOC activity
+figure(89)
+if strcmp(saveAN_spikesOrProbability,'probability')
+    % brainstem activity
+    PSTH=UTIL_PSTHmaker...
+        (ANprobRateOutput(2,:), dt, PSTHbinwidth)*dt/PSTHbinwidth;
+else
+    % AN probability
+    PSTH=UTIL_PSTHmaker(ICoutput(2,:), dtSpikes, PSTHbinwidth)*dt/PSTHbinwidth;
+end
+
+time=PSTHbinwidth:PSTHbinwidth:PSTHbinwidth*length(PSTH);
+plot(time, PSTH)
+set(gcf,'name', 'Lieberman')
+title(saveAN_spikesOrProbability)
+
+toc
+path(restorePath)
+
+% figure(88), plot(MOCattenuation)