rmeddis@29: function testOME(paramsName, paramChanges)
rmeddis@38: % testOME compute the external resonance and 
rmeddis@38: %  stapes response at a number of frequencies
rmeddis@38: % It compares the stapes displacement against human in vivo data
rmeddis@35: %  collected by Huber et al.2001.
rmeddis@38: % paramsName: name of file in parameterStore containing model parameters
rmeddis@38: % paramchanges: string array of changes to parameters. 
rmeddis@38: %    this can be omitted or {} is acceptable argument
rmeddis@38: % 
rmeddis@35: % testOME('Normal',{})
rmeddis@29: 
rmeddis@29: savePath=path;
rmeddis@29: addpath (['..' filesep 'utilities'],['..' filesep 'MAP'])
rmeddis@29: 
rmeddis@38: % default arguments
rmeddis@29: if nargin<2
rmeddis@29:     paramChanges=[];
rmeddis@29: end
rmeddis@38: if nargin<1
rmeddis@38:     paramsName='Normal';
rmeddis@38: end
rmeddis@29: 
rmeddis@29: sampleRate=50000;
rmeddis@29: 
rmeddis@29: dt=1/sampleRate;
rmeddis@29: leveldBSPL=80; 		% dB SPL as used by Huber (may trigger AR)
rmeddis@29: amp=10^(leveldBSPL/20)*28e-6;
rmeddis@29: duration=.05;
rmeddis@29: time=dt: dt: duration;
rmeddis@29: 
rmeddis@29: %% Comparison data (human)
rmeddis@29: % These data are taken directly from Huber 2001 (Fig. 4)
rmeddis@29: HuberFrequencies=[600	  800	 1000	 2000	 3000	4000 6000 8000];
rmeddis@29: HuberDisplacementAt80dBSPL=[1.5E-9	1.5E-09	1.5E-09	1.0E-09	7.0E-10	...
rmeddis@29:     3.0E-10	2.0E-10	1.0E-10]; % m;
rmeddis@29: % HuberVelocityAt80dBSPL= 2*pi*HuberFrequencies.*HuberDisplacementAt80dBSPL;
rmeddis@29: 
rmeddis@35: figure(2), clf, subplot(2,1,2)
rmeddis@29: set(2,'position',[5   349   268   327])
rmeddis@29: semilogx(HuberFrequencies, 20*log10(HuberDisplacementAt80dBSPL/1e-10),...
rmeddis@29:     'ko', 'MarkerFaceColor','k', 'Marker','o', 'markerSize',6)
rmeddis@29: hold on
rmeddis@29: 
rmeddis@29: %% Generate test stimulus .................................................................
rmeddis@29: 
rmeddis@29: % independent test using discrete frequencies
rmeddis@29: peakResponses=[];
rmeddis@29: peakTMpressure=[];
rmeddis@29: frequencies=[200 400 HuberFrequencies 10000];
rmeddis@29: for toneFrequency=frequencies
rmeddis@29:     inputSignal=amp*sin(2*pi*toneFrequency*time);
rmeddis@29: 
rmeddis@29:     showPlotsAndDetails=0;
rmeddis@29:     AN_spikesOrProbability='probability';
rmeddis@35:     
rmeddis@29:     % switch off AR & MOC (Huber's patients were deaf)
rmeddis@35:     idx=length(paramChanges);
rmeddis@35:     paramChanges{idx+1}='OMEParams.rateToAttenuationFactorProb=0;';
rmeddis@35:     paramChanges{idx+2}='DRNLParams.rateToAttenuationFactorProb = 0;';
rmeddis@29: 
rmeddis@29:     global OMEoutput  OMEextEarPressure TMoutput ARattenuation
rmeddis@29:     % BF is irrelevant
rmeddis@29:     MAP1_14(inputSignal, sampleRate, -1, ...
rmeddis@29:         paramsName, AN_spikesOrProbability, paramChanges);
rmeddis@29: 
rmeddis@29:     peakDisplacement=max(OMEoutput(floor(end/2):end));
rmeddis@29:     peakResponses=[peakResponses peakDisplacement];
rmeddis@29: 
rmeddis@29:     peakTMpressure=[peakTMpressure max(OMEextEarPressure)];
rmeddis@29: end
rmeddis@29: 
rmeddis@29: %% Report
rmeddis@29: disp('frequency displacement(m)')
rmeddis@29: % disp(num2str([frequencies' peakResponses']))
rmeddis@29: fprintf('%6.0f \t%10.3e\n',[frequencies' peakResponses']')
rmeddis@29: 
rmeddis@29: % stapes peak displacement
rmeddis@35: figure(2), subplot(2,1,2), hold on
rmeddis@29: semilogx(frequencies, 20*log10(peakResponses/1e-10), 'r', 'linewidth',4)
rmeddis@29: set(gca,'xScale','log')
rmeddis@29: % ylim([1e-11 1e-8])
rmeddis@29: xlim([100 10000]), ylim([0 30])
rmeddis@29: grid on
rmeddis@35: title(['stapes at ' num2str(leveldBSPL)])
rmeddis@29: ylabel('disp: dB re 1e-10m')
rmeddis@29: xlabel('stimulus frequency (Hz)')
rmeddis@29: legend({'Huber et al','model'},'location','southWest')
rmeddis@29: set(gcf,'name','OME')
rmeddis@29: 
rmeddis@29: % external ear resonance
rmeddis@35: figure(2), subplot(2,1,1),hold off
rmeddis@29: semilogx(frequencies, 20*log10(peakTMpressure/28e-6)-leveldBSPL,...
rmeddis@29:     'k', 'linewidth',2)
rmeddis@29: xlim([100 10000]) %, ylim([-10 30])
rmeddis@29: grid on
rmeddis@29: title(['External ear resonances' ])
rmeddis@29: ylabel('dB')
rmeddis@29: xlabel('frequency')
rmeddis@29: set(gcf,'name','OME: external resonances')
rmeddis@29: % ---------------------------------------------------------- display parameters
rmeddis@29: disp(['parameter file was: ' paramsName])
rmeddis@29: fprintf('\n')
rmeddis@29: 
rmeddis@29: path(savePath);