comparison testPrograms/test_DRNL_Ruggero97.m @ 38:c2204b18f4a2 tip

End nov big change

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