rmeddis@0: function showMAP (options) rmeddis@0: % defaults rmeddis@0: % options.showModelParameters=1; rmeddis@0: % options.showModelOutput=1; rmeddis@0: % options.printFiringRates=1; rmeddis@0: % options.showACF=1; rmeddis@0: % options.showEfferent=1; rmeddis@0: rmeddis@0: dbstop if warning rmeddis@0: rmeddis@0: global dt ANdt saveAN_spikesOrProbability savedBFlist saveMAPparamsName... rmeddis@0: savedInputSignal TMoutput OMEoutput ARattenuation ... rmeddis@0: DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV... rmeddis@0: IHCoutput ANprobRateOutput ANoutput savePavailable tauCas ... rmeddis@0: CNoutput ICoutput ICmembraneOutput ICfiberTypeRates MOCattenuation rmeddis@0: rmeddis@0: restorePath=path; rmeddis@0: addpath ( ['..' filesep 'utilities'], ['..' filesep 'parameterStore']) rmeddis@0: rmeddis@0: if nargin<1 rmeddis@0: options.showModelParameters=1; rmeddis@0: options.showModelOutput=1; rmeddis@0: options.printFiringRates=1; rmeddis@0: options.showACF=1; rmeddis@0: options.showEfferent=1; rmeddis@0: end rmeddis@0: rmeddis@0: if options.showModelParameters rmeddis@0: % Read parameters from MAPparams<***> file in 'parameterStore' folder rmeddis@0: % and print out all parameters rmeddis@0: cmd=['MAPparams' saveMAPparamsName ... rmeddis@0: '(-1, 1/dt, 1);']; rmeddis@0: eval(cmd); rmeddis@0: end rmeddis@0: rmeddis@0: if options.printFiringRates rmeddis@0: %% print summary firing rates rmeddis@0: fprintf('\n\n') rmeddis@0: disp('summary') rmeddis@0: disp(['AR: ' num2str(min(ARattenuation))]) rmeddis@0: disp(['MOC: ' num2str(min(min(MOCattenuation)))]) rmeddis@0: nANfiberTypes=length(tauCas); rmeddis@0: if strcmp(saveAN_spikesOrProbability, 'spikes') rmeddis@0: nANfibers=size(ANoutput,1); rmeddis@0: nHSRfibers=nANfibers/nANfiberTypes; rmeddis@0: duration=size(TMoutput,2)*dt; rmeddis@0: disp(['AN: ' num2str(sum(sum(ANoutput(end-nHSRfibers+1:end,:)))/... rmeddis@0: (nHSRfibers*duration))]) rmeddis@0: rmeddis@0: nCNneurons=size(CNoutput,1); rmeddis@0: nHSRCNneuronss=nCNneurons/nANfiberTypes; rmeddis@0: disp(['CN: ' num2str(sum(sum(CNoutput(end-nHSRCNneuronss+1:end,:)))... rmeddis@0: /(nHSRCNneuronss*duration))]) rmeddis@0: disp(['IC: ' num2str(sum(sum(ICoutput)))]) rmeddis@0: disp(['IC by type: ' num2str(mean(ICfiberTypeRates,2)')]) rmeddis@0: else rmeddis@0: disp(['AN: ' num2str(mean(mean(ANprobRateOutput)))]) rmeddis@0: end rmeddis@0: end rmeddis@0: rmeddis@0: rmeddis@0: %% figure (99) summarises main model output rmeddis@0: if options.showModelOutput rmeddis@0: plotInstructions.figureNo=99; rmeddis@0: signalRMS=mean(savedInputSignal.^2)^0.5; rmeddis@0: signalRMSdb=20*log10(signalRMS/20e-6); rmeddis@0: rmeddis@0: % plot signal (1) rmeddis@0: plotInstructions.displaydt=dt; rmeddis@0: plotInstructions.numPlots=6; rmeddis@0: plotInstructions.subPlotNo=1; rmeddis@0: plotInstructions.title=... rmeddis@0: ['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']; rmeddis@0: r=size(savedInputSignal,1); rmeddis@0: if r==1, savedInputSignal=savedInputSignal'; end rmeddis@0: UTIL_plotMatrix(savedInputSignal', plotInstructions); rmeddis@0: rmeddis@0: % stapes (2) rmeddis@0: plotInstructions.subPlotNo=2; rmeddis@0: plotInstructions.title= ['stapes displacement']; rmeddis@0: UTIL_plotMatrix(OMEoutput, plotInstructions); rmeddis@0: rmeddis@0: % DRNL (3) rmeddis@0: plotInstructions.subPlotNo=3; rmeddis@0: plotInstructions.title= ['BM displacement']; rmeddis@0: plotInstructions.yValues= savedBFlist; rmeddis@0: UTIL_plotMatrix(DRNLoutput, plotInstructions); rmeddis@0: rmeddis@0: switch saveAN_spikesOrProbability rmeddis@0: case 'spikes' rmeddis@0: % AN (4) rmeddis@0: plotInstructions.displaydt=ANdt; rmeddis@0: plotInstructions.title='AN'; rmeddis@0: plotInstructions.subPlotNo=4; rmeddis@0: plotInstructions.yLabel='BF'; rmeddis@0: plotInstructions.yValues= savedBFlist; rmeddis@0: plotInstructions.rasterDotSize=1; rmeddis@0: plotInstructions.plotDivider=1; rmeddis@0: if sum(sum(ANoutput))<100 rmeddis@0: plotInstructions.rasterDotSize=3; rmeddis@0: end rmeddis@0: UTIL_plotMatrix(ANoutput, plotInstructions); rmeddis@0: rmeddis@0: % CN (5) rmeddis@0: plotInstructions.displaydt=ANdt; rmeddis@0: plotInstructions.subPlotNo=5; rmeddis@0: plotInstructions.title='CN spikes'; rmeddis@0: if sum(sum(CNoutput))<100 rmeddis@0: plotInstructions.rasterDotSize=3; rmeddis@0: end rmeddis@0: UTIL_plotMatrix(CNoutput, plotInstructions); rmeddis@0: rmeddis@0: % IC (6) rmeddis@0: plotInstructions.displaydt=ANdt; rmeddis@0: plotInstructions.subPlotNo=6; rmeddis@0: plotInstructions.title='IC'; rmeddis@0: if size(ICoutput,1)>3 rmeddis@0: if sum(sum(ICoutput))<100 rmeddis@0: plotInstructions.rasterDotSize=3; rmeddis@0: end rmeddis@0: UTIL_plotMatrix(ICoutput, plotInstructions); rmeddis@0: else rmeddis@0: plotInstructions.title='IC (HSR) membrane potential'; rmeddis@0: plotInstructions.displaydt=dt; rmeddis@0: plotInstructions.yLabel='V'; rmeddis@0: plotInstructions.zValuesRange= [-.1 0]; rmeddis@0: UTIL_plotMatrix(ICmembraneOutput, plotInstructions); rmeddis@0: end rmeddis@0: rmeddis@0: otherwise % probability (4-6) rmeddis@0: plotInstructions.displaydt=dt; rmeddis@0: plotInstructions.numPlots=2; rmeddis@0: plotInstructions.subPlotNo=2; rmeddis@0: plotInstructions.yLabel='BF'; rmeddis@0: if nANfiberTypes>1, rmeddis@0: plotInstructions.yLabel='LSR HSR'; rmeddis@0: plotInstructions.plotDivider=1; rmeddis@0: end rmeddis@0: plotInstructions.title='AN - spike probability'; rmeddis@0: UTIL_plotMatrix(ANprobRateOutput, plotInstructions); rmeddis@0: end rmeddis@0: end rmeddis@0: rmeddis@0: %% plot efferent control values as dB rmeddis@0: if options.showEfferent rmeddis@0: plotInstructions=[]; rmeddis@0: plotInstructions.figureNo=98; rmeddis@0: figure(98), clf rmeddis@0: plotInstructions.displaydt=dt; rmeddis@0: plotInstructions.numPlots=2; rmeddis@0: plotInstructions.subPlotNo=1; rmeddis@0: plotInstructions.zValuesRange=[ -25 0]; rmeddis@0: plotInstructions.title= ['AR strength. Signal level= ' ... rmeddis@0: num2str(signalRMSdb,'%4.0f') ' dB SPL']; rmeddis@0: UTIL_plotMatrix(20*log10(ARattenuation), plotInstructions); rmeddis@0: rmeddis@0: plotInstructions.subPlotNo=2; rmeddis@0: plotInstructions.yValues= savedBFlist; rmeddis@0: plotInstructions.yLabel= 'BF'; rmeddis@0: plotInstructions.title= ['MOC strength']; rmeddis@0: plotInstructions.zValuesRange=[ -25 0]; rmeddis@0: subplot(2,1,2) rmeddis@0: % imagesc(MOCattenuation) rmeddis@0: UTIL_plotMatrix(20*log10(MOCattenuation), plotInstructions); rmeddis@0: colorbar rmeddis@0: end rmeddis@0: rmeddis@0: rmeddis@0: %% ACF plot if required rmeddis@0: if options.showACF rmeddis@0: tic rmeddis@0: method.dt=dt; rmeddis@0: method.segmentNo=1; rmeddis@0: method.nonlinCF=savedBFlist; rmeddis@0: rmeddis@0: minPitch= 80; maxPitch= 4000; numPitches=100; % specify lags rmeddis@0: pitches=10.^ linspace(log10(minPitch), log10(maxPitch),numPitches); rmeddis@0: pitches=fliplr(pitches); rmeddis@0: filteredSACFParams.lags=1./pitches; % autocorrelation lags vector rmeddis@0: filteredSACFParams.acfTau= .003; % time constant of running ACF rmeddis@0: filteredSACFParams.lambda= 0.12; % slower filter to smooth ACF rmeddis@0: filteredSACFParams.lambda= 0.01; % slower filter to smooth ACF rmeddis@0: rmeddis@0: filteredSACFParams.plotACFs=0; % special plot (see code) rmeddis@0: filteredSACFParams.plotFilteredSACF=0; % 0 plots unfiltered ACFs rmeddis@0: filteredSACFParams.plotMoviePauses=.3; % special plot (see code) rmeddis@0: rmeddis@0: filteredSACFParams.usePressnitzer=0; % attenuates ACF at long lags rmeddis@0: filteredSACFParams.lagsProcedure= 'useAllLags'; rmeddis@0: % filteredSACFParams.lagsProcedure= 'useBernsteinLagWeights'; rmeddis@0: % filteredSACFParams.lagsProcedure= 'omitShortLags'; rmeddis@0: filteredSACFParams.criterionForOmittingLags=3; rmeddis@0: filteredSACFParams.plotACFsInterval=200; rmeddis@0: rmeddis@0: if filteredSACFParams.plotACFs rmeddis@0: % plot original waveform on ACF plot rmeddis@0: figure(13), clf rmeddis@0: subplot(4,1,1) rmeddis@0: t=dt*(1:length(savedInputSignal)); rmeddis@0: plot(t,savedInputSignal) rmeddis@0: xlim([0 t(end)]) rmeddis@0: title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']); rmeddis@0: end rmeddis@0: rmeddis@0: % plot original waveform on summary/smoothed ACF plot rmeddis@0: figure(97), clf rmeddis@0: subplot(2,1,1) rmeddis@0: t=dt*(1:length(savedInputSignal)); rmeddis@0: plot(t,savedInputSignal) rmeddis@0: xlim([0 t(end)]) rmeddis@0: title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']); rmeddis@0: rmeddis@0: rmeddis@0: % compute ACF rmeddis@0: switch saveAN_spikesOrProbability rmeddis@0: case 'probability' rmeddis@0: inputToACF=ANprobRateOutput.^0.5; rmeddis@0: otherwise rmeddis@0: inputToACF=ANoutput; rmeddis@0: end rmeddis@0: rmeddis@0: disp ('computing ACF...') rmeddis@0: [P, BFlist, sacf, boundaryValue] = ... rmeddis@0: filteredSACF(inputToACF, method, filteredSACFParams); rmeddis@0: disp(' ACF done.') rmeddis@0: rmeddis@0: % SACF rmeddis@0: subplot(2,1,2) rmeddis@0: imagesc(P) rmeddis@0: ylabel('periodicities (Hz)') rmeddis@0: xlabel('time (s)') rmeddis@0: title(['running smoothed (root) SACF. ' saveAN_spikesOrProbability ' input']) rmeddis@0: pt=[1 get(gca,'ytick')]; % force top xtick to show rmeddis@0: set(gca,'ytick',pt) rmeddis@0: set(gca,'ytickLabel', round(pitches(pt))) rmeddis@0: tt=get(gca,'xtick'); rmeddis@0: set(gca,'xtickLabel', round(100*t(tt))/100) rmeddis@0: end rmeddis@0: rmeddis@0: path(restorePath)