rmeddis@38: function test_DRNL_Ruggero97
rmeddis@38: % test_DRNL_Ruggero97 attempts to match Ruggero's (1992 and 1997)
rmeddis@38: %  iso-intensity data by fiddling with the parameters
rmeddis@38: 
rmeddis@38: % # BF is the BF of the filter to be assessed
rmeddis@38: BF=9000;
rmeddis@38: 
rmeddis@38: % # test frequencies. check that BF is one of them
rmeddis@38: %    copy Ruggero's test tones as fara as possible
rmeddis@38: numFs=6; lowestF=4000; 	highestF= 11030;
rmeddis@38: toneFrequencyList=round(logspace(log10(lowestF), log10(highestF), numFs));
rmeddis@38: 
rmeddis@38: %  # parameter file name. this is the base set of parameters
rmeddis@38: MAPparamsName='Normal';
rmeddis@38: 
rmeddis@38: % # probability representation (not directly relevant here as only
rmeddis@38: %    the DRNL output is used
rmeddis@38: AN_spikesOrProbability='probability';
rmeddis@38: 
rmeddis@38: % # tone characteristics
rmeddis@38: sampleRate= 100000;
rmeddis@38: duration=0.0200;                % Ruggero uses 5, 10, 25 ms tones
rmeddis@38: rampDuration=0.0015;              % raised cosine ramp (seconds)
rmeddis@38: beginSilence=0.050;
rmeddis@38: endSilence=0.020;
rmeddis@38: 
rmeddis@38: % # levels
rmeddis@38: levels=[3 10:10:80];
rmeddis@38: 
rmeddis@38: %% # change model parameters
rmeddis@38: % Parameter changes can be used to change one or more model parameters
rmeddis@38: %  *after* the MAPparams file has been read
rmeddis@38: 
rmeddis@38: paramChanges={};
rmeddis@38: % switch off all efferent effects and then play with DRNL params
rmeddis@38: paramChanges={...
rmeddis@38:     'DRNLParams.rateToAttenuationFactorProb = 0.00; ',...
rmeddis@38:     'OMEParams.rateToAttenuationFactorProb=0.0;', ...
rmeddis@38:     'DRNLParams.ctBMdB = -20;'...
rmeddis@38:     'DRNLParams.g=1000;'...
rmeddis@38:     'DRNLParams.linCFs=7000;'...
rmeddis@38:     'DRNLParams.linBWs=3500;'...
rmeddis@38:     };
rmeddis@38: 
rmeddis@38: % delare 'showMap' options to control graphical output
rmeddis@38: % (not needed but might be useful
rmeddis@38: % showMapOptions.printModelParameters=0;   % prints all parameters
rmeddis@38: % showMapOptions.showModelOutput=1;       % plot of all stages
rmeddis@38: % showMapOptions.printFiringRates=0;      % prints stage activity levels
rmeddis@38: % showMapOptions.showACF=0;               % shows SACF (probability only)
rmeddis@38: % showMapOptions.showEfferent=0;          % tracks of AR and MOC
rmeddis@38: % showMapOptions.surfProbability=0;       % 2D plot of HSR response
rmeddis@38: % showMapOptions.surfSpikes=0;            % 2D plot of spikes histogram
rmeddis@38: % showMapOptions.ICrates=0;               % IC rates by CNtauGk
rmeddis@38: % showMapOptions.fileName=[];
rmeddis@38: 
rmeddis@38: % run the program
rmeddis@38: global dt  DRNLoutput DRNLParams TMoutput
rmeddis@38: dbstop if error
rmeddis@38: restorePath=path; 
rmeddis@38: addpath (['..' filesep 'MAP'],    ['..' filesep 'wavFileStore'], ...
rmeddis@38:     ['..' filesep 'utilities'])
rmeddis@38: figure(2), clf,figure(3), clf,figure(4), clf,
rmeddis@38: 
rmeddis@38: peakOutputDisp=zeros(length(levels),length(toneFrequencyList));
rmeddis@38: peakStapesDisp=zeros(length(levels),length(toneFrequencyList));
rmeddis@38: 
rmeddis@38: %% now vary level and frequency while measuring the response
rmeddis@38: levelNo=0;
rmeddis@38: for leveldBSPL=levels
rmeddis@38:     levelNo=levelNo+1;
rmeddis@38:     disp(['level: ' num2str(leveldBSPL)])
rmeddis@38: 
rmeddis@38:     freqNo=0;
rmeddis@38:     for toneFrequency=toneFrequencyList
rmeddis@38:         freqNo=freqNo+1;
rmeddis@38: 
rmeddis@38:         %% Generate stimuli
rmeddis@38:         dt=1/sampleRate;
rmeddis@38:         time=dt: dt: duration;
rmeddis@38:         inputSignal=sum(sin(2*pi*toneFrequency'*time), 1);
rmeddis@38:         amp=10^(leveldBSPL/20)*28e-6;   % converts to Pascals (peak)
rmeddis@38:         inputSignal=amp*inputSignal;
rmeddis@38:         % apply ramps
rmeddis@38:         if rampDuration>0.5*duration, rampDuration=duration/2; end
rmeddis@38:         rampTime=dt:dt:rampDuration;
rmeddis@38:         ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
rmeddis@38:             ones(1,length(time)-length(rampTime))];
rmeddis@38:         inputSignal=inputSignal.*ramp;
rmeddis@38:         ramp=fliplr(ramp);
rmeddis@38:         inputSignal=inputSignal.*ramp;
rmeddis@38:         % add silence
rmeddis@38:         intialSilence= zeros(1,round(beginSilence/dt));
rmeddis@38:         finalSilence= zeros(1,round(endSilence/dt));
rmeddis@38:         inputSignal= [intialSilence inputSignal finalSilence];
rmeddis@38: 
rmeddis@38:         %% run the model
rmeddis@38:     
rmeddis@38:         MAP1_14(inputSignal, sampleRate, BF, ...
rmeddis@38:             MAPparamsName, AN_spikesOrProbability, paramChanges);
rmeddis@38: 
rmeddis@38:         peakOutputDisp(levelNo,freqNo)=max(DRNLoutput);
rmeddis@38:         peakStapesDisp(levelNo,freqNo)=max(TMoutput);
rmeddis@38: 
rmeddis@38:         % the model run is complete. Now display the results if debugging
rmeddis@38:         % UTIL_showMAP(showMapOptions, paramChanges)
rmeddis@38:     end % probe frequencies
rmeddis@38:     
rmeddis@38:     % monitor progress
rmeddis@38:     figure(2), loglog(toneFrequencyList, peakOutputDisp), hold on
rmeddis@38:     xlabel('frequency')
rmeddis@38:     ylabel('peak DRNL displacement (m)')
rmeddis@38:     title ('DRNL uncorrected displacement')
rmeddis@38: end  % levels
rmeddis@38: figure(2),legend(num2str(toneFrequencyList'),'location','northwest')
rmeddis@38: 
rmeddis@38: % convert from model BM displacement to Ruggero's velocity
rmeddis@38: DRNLvelocity= peakOutputDisp ...
rmeddis@38:     .*repmat(2*pi*toneFrequencyList,length(levels),1);
rmeddis@38: % convert from model stapes displacement to Ruggero's velocity
rmeddis@38: stapesVelocity= peakStapesDisp ...
rmeddis@38:     .*repmat(2*pi*toneFrequencyList,length(levels),1);
rmeddis@38: % velocity gain is increased velocity attributable to the DRNL
rmeddis@38: DRNLvelocityGain = DRNLvelocity./stapesVelocity;
rmeddis@38: DRNLvelocityGain_dB= 20*log10(DRNLvelocityGain );
rmeddis@38: % iso-intensity equates all functions at the same outpu for the lowest
rmeddis@38: %  velocity tested
rmeddis@38: isoIntensityDRNLvel_dB= DRNLvelocityGain_dB- ...
rmeddis@38:     repmat(DRNLvelocityGain_dB(:,1),1,length(toneFrequencyList));
rmeddis@38: 
rmeddis@38: % displays
rmeddis@38: % iso-velocity function (nrmalised by stapes velocity)
rmeddis@38: figure(3), clf, semilogx(toneFrequencyList, isoIntensityDRNLvel_dB)
rmeddis@38: ylim([-10 60])
rmeddis@38: xlim([2000 20000])
rmeddis@38: xlabel('frequency')
rmeddis@38: ylabel('peak DRNL velocity gain (dB)')
rmeddis@38: title(['CF= ' num2str(BF) ' Hz'])
rmeddis@38: legend(num2str(levels'),'location','northwest')
rmeddis@38: 
rmeddis@38: % velocity I/O function
rmeddis@38: figure(4), clf, semilogy(levels, (DRNLvelocity*1e6)'), hold on
rmeddis@38: ylim([1e0 1e4])
rmeddis@38: xlim([0 100])
rmeddis@38: xlabel('level')
rmeddis@38: ylabel(' velocity (microm/ s)')
rmeddis@38: title(['CF= ' num2str(BF) ' Hz'])
rmeddis@38: legend(num2str(toneFrequencyList'),'location','northwest')
rmeddis@38: 
rmeddis@38: % command window reports
rmeddis@38: disp(''), disp('***')
rmeddis@38: disp(AN_spikesOrProbability)
rmeddis@38: 
rmeddis@38: % DRNL parameter set
rmeddis@38: UTIL_showStructureSummary(DRNLParams, 'DRNLParams', 10)
rmeddis@38: 
rmeddis@38: % stimulus characteristics
rmeddis@38: disp(['CF=' num2str(BF)])
rmeddis@38: disp(['tone Duration=' num2str(rampDuration)])
rmeddis@38: disp(['ramp Duration=' num2str(duration)])
rmeddis@38: 
rmeddis@38: % particular parameter changes used on this run
rmeddis@38: for i=1:length(paramChanges)
rmeddis@38:     disp(paramChanges{i})
rmeddis@38: end
rmeddis@38: 
rmeddis@38: % if required dump one or more matrices in tab spaced format
rmeddis@38: %  (suitable for pasting directly into EXCEL.
rmeddis@38: %     UTIL_printTabTable([toneFrequencyList' peakOutputDisp'])
rmeddis@38: 
rmeddis@38: % leave the path as you found it
rmeddis@38: path(restorePath)
rmeddis@38: 
rmeddis@38: