annotate 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)