view testPrograms/testOME.m @ 38:c2204b18f4a2 tip

End nov big change

function testOME(paramsName, paramChanges)
% testOME compute the external resonance and 
%  stapes response at a number of frequencies
% It compares the stapes displacement against human in vivo data
%  collected by Huber et al.2001.
% paramsName: name of file in parameterStore containing model parameters
% paramchanges: string array of changes to parameters. 
%    this can be omitted or {} is acceptable argument
% 
% testOME('Normal',{})

savePath=path;
addpath (['..' filesep 'utilities'],['..' filesep 'MAP'])

% default arguments
if nargin<2
    paramChanges=[];
end
if nargin<1
    paramsName='Normal';
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,2)
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)
    idx=length(paramChanges);
    paramChanges{idx+1}='OMEParams.rateToAttenuationFactorProb=0;';
    paramChanges{idx+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)];
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,2), 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)])
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,1),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);