Mercurial > hg > map
changeset 19:5b23b9f11806
MAP debugging
| author | Ray Meddis <rmeddis@essex.ac.uk> |
|---|---|
| date | Tue, 07 Jun 2011 17:50:05 +0100 |
| parents | e9e263e4fcde |
| children | fafe69c43108 |
| files | MAP/MAP1_14.m multithreshold 1.46/savedData/mostRecentResults.mat multithreshold 1.46/testBM.m multithreshold 1.46/testOME.m parameterStore/MAPparamsNormal.m testPrograms/showMAP.m testPrograms/test_MAP1_14.m utilities/UTIL_makePSTH.m |
| diffstat | 8 files changed, 179 insertions(+), 142 deletions(-) [+] |
line wrap: on
line diff
--- a/MAP/MAP1_14.m Tue Jun 07 10:54:42 2011 +0100 +++ b/MAP/MAP1_14.m Tue Jun 07 17:50:05 2011 +0100 @@ -535,34 +535,34 @@ % shorter segments after speed up. shorterSegmentEndPTR=reducedSegmentPTR+reducedSegmentLength-1; - iputPressureSegment=inputSignal... + inputPressureSegment=inputSignal... (:,segmentStartPTR:segmentStartPTR+segmentLength-1); % segment debugging plots % figure(98) - % plot(segmentTime,iputPressureSegment), title('signalSegment') + % plot(segmentTime,inputPressureSegment), title('signalSegment') % OME ---------------------- % OME Stage 1: external resonances. Add to inputSignal pressure wave - y=iputPressureSegment; + y=inputPressureSegment; for n=1:nOMEExtFilters % any number of resonances can be used [x OMEExtFilterBndry{n}] = ... filter(ExtFilter_b{n},ExtFilter_a{n},... - iputPressureSegment, OMEExtFilterBndry{n}); + inputPressureSegment, OMEExtFilterBndry{n}); x= x* OMEgainScalars(n); % This is a parallel resonance so add it y=y+x; end - iputPressureSegment=y; - OMEextEarPressure(segmentStartPTR:segmentEndPTR)= iputPressureSegment; + inputPressureSegment=y; + OMEextEarPressure(segmentStartPTR:segmentEndPTR)= inputPressureSegment; % OME stage 2: convert input pressure (velocity) to % tympanic membrane(TM) displacement using low pass filter [TMdisplacementSegment OME_TMdisplacementBndry] = ... - filter(TMdisp_b,TMdisp_a,iputPressureSegment, ... + filter(TMdisp_b,TMdisp_a,inputPressureSegment, ... OME_TMdisplacementBndry); % and save it TMoutput(segmentStartPTR:segmentEndPTR)= TMdisplacementSegment; @@ -610,6 +610,7 @@ else % no MOC available yet MOC=ones(1, segmentLength); end + plot(MOC) % current channel % first gammatone filter for order = 1 : GTnonlinOrder @@ -623,10 +624,10 @@ y= nonlinOutput.*(MOC* DRNLa); % linear section. % compress those parts of the signal above the compression % threshold - abs_x = abs(nonlinOutput); + abs_x = abs(y); idx=find(abs_x>DRNLcompressionThreshold); if ~isempty(idx)>0 - y(idx)=sign(nonlinOutput(idx)).*... + y(idx)=sign(y(idx)).*... (DRNLb*abs_x(idx).^DRNLc); end nonlinOutput=y; @@ -634,7 +635,8 @@ % second filter removes distortion products for order = 1 : GTnonlinOrder [ nonlinOutput GTnonlinBdry2{BFno,order}] = ... - filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), nonlinOutput, GTnonlinBdry2{BFno,order}); + filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ... + nonlinOutput, GTnonlinBdry2{BFno,order}); end % combine the two paths to give the DRNL displacement @@ -757,14 +759,18 @@ % Estimate efferent effects. ARattenuation (0 <> 1) % acoustic reflex - ARAttSeg=mean(ANrate(1:nBFs,:),1); %LSR channels are 1:nBF - % smooth - [ARAttSeg, ARboundaryProb] = ... - filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundaryProb); - ARAttSeg=ARAttSeg-ARrateThreshold; - ARAttSeg(ARAttSeg<0)=0; % prevent negative strengths - ARattenuation(segmentStartPTR:segmentEndPTR)=... - (1-ARrateToAttenuationFactorProb.* ARAttSeg); + [r c]=size(ANrate); + if r>nBFs % Only if LSR fibers are computed + ARAttSeg=mean(ANrate(1:nBFs,:),1); %LSR channels are 1:nBF + % smooth + [ARAttSeg, ARboundaryProb] = ... + filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundaryProb); + ARAttSeg=ARAttSeg-ARrateThreshold; + ARAttSeg(ARAttSeg<0)=0; % prevent negative strengths + ARattenuation(segmentStartPTR:segmentEndPTR)=... + (1-ARrateToAttenuationFactorProb.* ARAttSeg); + end + % plot(ARattenuation) % MOC attenuation % within-channel HSR response only @@ -780,8 +786,8 @@ (1- smoothedRates* rateToAttenuationFactorProb); end MOCattenuation(MOCattenuation<0)=0.001; -plot(MOCattenuation) -MOCattenuation + + % plot(MOCattenuation) case 'spikes'
--- a/multithreshold 1.46/testBM.m Tue Jun 07 10:54:42 2011 +0100 +++ b/multithreshold 1.46/testBM.m Tue Jun 07 17:50:05 2011 +0100 @@ -19,7 +19,7 @@ % levels= 50; nLevels=length(levels); relativeFrequencies=[0.25 .5 .75 1 1.25 1.5 2]; -% relativeFrequencies=1; +relativeFrequencies=1; % refBMdisplacement is the displacement of the BM at threshold % 1 nm disp at threshold (9 kHz, Ruggero)
--- a/multithreshold 1.46/testOME.m Tue Jun 07 10:54:42 2011 +0100 +++ b/multithreshold 1.46/testOME.m Tue Jun 07 17:50:05 2011 +0100 @@ -23,9 +23,11 @@ set(2,'position',[5 349 268 327]) semilogx(HuberFrequencies, 20*log10(HuberDisplacementAt80dBSPL/1e-10),... 'ko', 'MarkerFaceColor','k', 'Marker','o', 'markerSize',6) -% Generate test stimulus ................................................................. +hold on -%% independent test using discrete frequencies +%% Generate test stimulus ................................................................. + +% independent test using discrete frequencies peakResponses=[]; peakTMpressure=[]; frequencies=[200 400 HuberFrequencies 10000]; @@ -39,7 +41,7 @@ paramChanges{1}='OMEParams.rateToAttenuationFactorProb=0;'; paramChanges{2}='DRNLParams.rateToAttenuationFactorProb = 0;'; - global OMEoutput OMEextEarPressure TMoutput ARAttenuation + global OMEoutput OMEextEarPressure TMoutput ARattenuation % BF is irrelevant MAP1_14(inputSignal, sampleRate, -1, ... MAPparamsName, AN_spikesOrProbability, paramChanges); @@ -48,7 +50,7 @@ peakResponses=[peakResponses peakDisplacement]; peakTMpressure=[peakTMpressure max(OMEextEarPressure)]; - disp([' greatest AR attenuation: ' num2str(min(ARAttenuation))]) + disp([' AR attenuation (dB): ' num2str(20*log10(min(ARattenuation)))]) end %% Report @@ -59,6 +61,7 @@ % stapes peak displacement figure(2), subplot(2,1,1), 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
--- a/parameterStore/MAPparamsNormal.m Tue Jun 07 10:54:42 2011 +0100 +++ b/parameterStore/MAPparamsNormal.m Tue Jun 07 17:50:05 2011 +0100 @@ -57,7 +57,7 @@ % Acoustic reflex: maximum attenuation should be around 25 dB Price (1966) % i.e. a minimum ratio of 0.056. % 'spikes' model: AR based on brainstem spiking activity (LSR) -OMEParams.rateToAttenuationFactor=0.003; % * N(all ICspikes) +OMEParams.rateToAttenuationFactor=0.004; % * N(all ICspikes) % OMEParams.rateToAttenuationFactor=0; % * N(all ICspikes) % 'probability model': Ar based on AN firing probabilities (LSR) @@ -76,7 +76,7 @@ % DRNL nonlinear path DRNLParams.a=3e4; % nonlinear path gain (below compression threshold) -% DRNLParams.a=3e2; % DRNL.a=0 means no OHCs (no nonlinear path) +DRNLParams.a=5e2; % DRNL.a=0 means no OHCs (no nonlinear path) DRNLParams.b=8e-6; % *compression threshold raised compression % DRNLParams.b=1; % b=1 means no compression @@ -102,13 +102,13 @@ DRNLParams.MOCdelay = efferentDelay; % must be < segment length! % 'spikes' model: MOC based on brainstem spiking activity (HSR) DRNLParams.rateToAttenuationFactor = .009; % strength of MOC -DRNLParams.rateToAttenuationFactor = .009; % strength of MOC +DRNLParams.rateToAttenuationFactor = .004; % strength of MOC % DRNLParams.rateToAttenuationFactor = 0; % strength of MOC % 'probability' model: MOC based on AN spiking activity (HSR) -DRNLParams.rateToAttenuationFactorProb = .007; % strength of MOC -DRNLParams.rateToAttenuationFactorProb = .0; % strength of MOC -DRNLParams.MOCtau =.03; % smoothing for MOC +DRNLParams.rateToAttenuationFactorProb = .004; % strength of MOC +% DRNLParams.rateToAttenuationFactorProb = .0; % strength of MOC +DRNLParams.MOCtau =.1; % smoothing for MOC DRNLParams.MOCrateThreshold =50; % set to AN rate threshold @@ -147,7 +147,7 @@ % reminder: changing z has a strong effect on HF thresholds (like Et) IHCpreSynapseParams.z= 2e42; % scalar Ca -> vesicle release rate -LSRtauCa=50e-6; HSRtauCa=85e-6; % seconds +LSRtauCa=35e-6; HSRtauCa=85e-6; % seconds % LSRtauCa=35e-6; HSRtauCa=70e-6; % seconds IHCpreSynapseParams.tauCa= [LSRtauCa HSRtauCa]; %LSR and HSR fiber
--- a/testPrograms/showMAP.m Tue Jun 07 10:54:42 2011 +0100 +++ b/testPrograms/showMAP.m Tue Jun 07 17:50:05 2011 +0100 @@ -53,7 +53,7 @@ duration=size(TMoutput,2)*dt; disp(['AN: ' num2str(sum(sum(ANoutput(end-nHSRfibers+1:end,:)))/... (nHSRfibers*duration))]) - + nCNneurons=size(CNoutput,1); nHSRCNneuronss=nCNneurons/nANfiberTypes; disp(['CN: ' num2str(sum(sum(CNoutput(end-nHSRCNneuronss+1:end,:)))... @@ -62,6 +62,9 @@ % disp(['IC by type: ' num2str(mean(ICfiberTypeRates,2)')]) else disp(['AN: ' num2str(mean(mean(ANprobRateOutput)))]) + [PSTH pointsPerBin]= UTIL_makePSTH(ANprobRateOutput, dt, 0.001); + disp(['max max AN: ' num2str(max(max(... + PSTH/pointsPerBin )))]) end end @@ -71,7 +74,7 @@ plotInstructions.figureNo=99; signalRMS=mean(savedInputSignal.^2)^0.5; signalRMSdb=20*log10(signalRMS/20e-6); - + % plot signal (1) plotInstructions.displaydt=dt; plotInstructions.numPlots=6; @@ -81,18 +84,18 @@ r=size(savedInputSignal,1); if r==1, savedInputSignal=savedInputSignal'; end UTIL_plotMatrix(savedInputSignal', plotInstructions); - + % stapes (2) plotInstructions.subPlotNo=2; plotInstructions.title= ['stapes displacement']; UTIL_plotMatrix(OMEoutput, plotInstructions); - + % DRNL (3) plotInstructions.subPlotNo=3; plotInstructions.title= ['BM displacement']; plotInstructions.yValues= savedBFlist; UTIL_plotMatrix(DRNLoutput, plotInstructions); - + switch saveAN_spikesOrProbability case 'spikes' % AN (4) @@ -107,7 +110,7 @@ plotInstructions.rasterDotSize=3; end UTIL_plotMatrix(ANoutput, plotInstructions); - + % CN (5) plotInstructions.displaydt=ANdt; plotInstructions.subPlotNo=5; @@ -116,7 +119,7 @@ plotInstructions.rasterDotSize=3; end UTIL_plotMatrix(CNoutput, plotInstructions); - + % IC (6) plotInstructions.displaydt=ANdt; plotInstructions.subPlotNo=6; @@ -133,7 +136,7 @@ plotInstructions.zValuesRange= [-.1 0]; UTIL_plotMatrix(ICmembraneOutput, plotInstructions); end - + otherwise % probability (4-6) plotInstructions.displaydt=dt; plotInstructions.numPlots=2; @@ -148,6 +151,32 @@ end end +%% surface plot of probability +if options.surfProbability + figure(97), clf + % select only HSR fibers at the bottom of the matrix + ANprobRateOutput= ANprobRateOutput(end-length(savedBFlist)+1:end,:); + [nY nX]=size(ANprobRateOutput); + if nY>2 + time=dt*(1:nX); + surf(time, savedBFlist, ANprobRateOutput) + shading interp + set(gca, 'yScale','log') + xlim([0 max(time)]), ylim([0 max(savedBFlist)]), zlim([0 1000]) + xlabel('time (s)') + ylabel('best frequency (Hz)') + zlabel('spike rate') + view([-20 60]) + if isfield(options, 'fileName') + title ([options.fileName ': ' num2str(signalRMSdb,'% 3.0f') ' dB']) + else + title ([ num2str(signalRMSdb,'% 3.0f') ' dB']) + end + + end +end + + %% plot efferent control values as dB if options.showEfferent plotInstructions=[]; @@ -160,7 +189,7 @@ plotInstructions.title= ['AR strength. Signal level= ' ... num2str(signalRMSdb,'%4.0f') ' dB SPL']; UTIL_plotMatrix(20*log10(ARattenuation), plotInstructions); - + plotInstructions.subPlotNo=2; plotInstructions.yValues= savedBFlist; plotInstructions.yLabel= 'BF'; @@ -171,95 +200,76 @@ UTIL_plotMatrix(20*log10(MOCattenuation), plotInstructions); colorbar end - -%% surface plot of probability -if options.surfProbability - figure(97), clf - % select only HSR fibers at the bottom of the matrix - ANprobRateOutput= ANprobRateOutput(end-length(savedBFlist)+1:end,:); - [nY nX]=size(ANprobRateOutput); - time=dt*(1:nX); - surf(time, savedBFlist, ANprobRateOutput) - shading interp - set(gca, 'yScale','log') - xlim([0 max(time)]), ylim([0 max(savedBFlist)]), zlim([0 1000]) - xlabel('time (s)') - ylabel('best frequency (Hz)') - zlabel('spike rate') - view([-20 60]) - title ([options.fileName ': ' num2str(signalRMSdb,'% 3.0f') ' dB']) -end - - -%% ACF plot if required -if options.showACF - tic - method.dt=dt; - method.segmentNo=1; - method.nonlinCF=savedBFlist; - - minPitch= 80; maxPitch= 4000; numPitches=100; % specify lags - pitches=10.^ linspace(log10(minPitch), log10(maxPitch),numPitches); - pitches=fliplr(pitches); - filteredSACFParams.lags=1./pitches; % autocorrelation lags vector - filteredSACFParams.acfTau= .003; % time constant of running ACF - filteredSACFParams.lambda= 0.12; % slower filter to smooth ACF - filteredSACFParams.lambda= 0.01; % slower filter to smooth ACF - - filteredSACFParams.plotACFs=0; % special plot (see code) - filteredSACFParams.plotFilteredSACF=0; % 0 plots unfiltered ACFs - filteredSACFParams.plotMoviePauses=.3; % special plot (see code) - - filteredSACFParams.usePressnitzer=0; % attenuates ACF at long lags - filteredSACFParams.lagsProcedure= 'useAllLags'; - % filteredSACFParams.lagsProcedure= 'useBernsteinLagWeights'; - % filteredSACFParams.lagsProcedure= 'omitShortLags'; - filteredSACFParams.criterionForOmittingLags=3; - filteredSACFParams.plotACFsInterval=200; - - if filteredSACFParams.plotACFs - % plot original waveform on ACF plot - figure(13), clf - subplot(4,1,1) + + %% ACF plot if required + if options.showACF + tic + method.dt=dt; + method.segmentNo=1; + method.nonlinCF=savedBFlist; + + minPitch= 80; maxPitch= 4000; numPitches=100; % specify lags + pitches=10.^ linspace(log10(minPitch), log10(maxPitch),numPitches); + pitches=fliplr(pitches); + filteredSACFParams.lags=1./pitches; % autocorrelation lags vector + filteredSACFParams.acfTau= .003; % time constant of running ACF + filteredSACFParams.lambda= 0.12; % slower filter to smooth ACF + filteredSACFParams.lambda= 0.01; % slower filter to smooth ACF + + filteredSACFParams.plotACFs=0; % special plot (see code) + filteredSACFParams.plotFilteredSACF=0; % 0 plots unfiltered ACFs + filteredSACFParams.plotMoviePauses=.3; % special plot (see code) + + filteredSACFParams.usePressnitzer=0; % attenuates ACF at long lags + filteredSACFParams.lagsProcedure= 'useAllLags'; + % filteredSACFParams.lagsProcedure= 'useBernsteinLagWeights'; + % filteredSACFParams.lagsProcedure= 'omitShortLags'; + filteredSACFParams.criterionForOmittingLags=3; + filteredSACFParams.plotACFsInterval=200; + + if filteredSACFParams.plotACFs + % plot original waveform on ACF plot + figure(13), clf + subplot(4,1,1) + t=dt*(1:length(savedInputSignal)); + plot(t,savedInputSignal) + xlim([0 t(end)]) + title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']); + end + + % plot original waveform on summary/smoothed ACF plot + figure(96), clf + subplot(2,1,1) t=dt*(1:length(savedInputSignal)); plot(t,savedInputSignal) xlim([0 t(end)]) title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']); + + + % compute ACF + switch saveAN_spikesOrProbability + case 'probability' + inputToACF=ANprobRateOutput.^0.5; + otherwise + inputToACF=ANoutput; + end + + disp ('computing ACF...') + [P, BFlist, sacf, boundaryValue] = ... + filteredSACF(inputToACF, method, filteredSACFParams); + disp(' ACF done.') + + % SACF + subplot(2,1,2) + imagesc(P) + ylabel('periodicities (Hz)') + xlabel('time (s)') + title(['running smoothed (root) SACF. ' saveAN_spikesOrProbability ' input']) + pt=[1 get(gca,'ytick')]; % force top xtick to show + set(gca,'ytick',pt) + set(gca,'ytickLabel', round(pitches(pt))) + tt=get(gca,'xtick'); + set(gca,'xtickLabel', round(100*t(tt))/100) end - - % plot original waveform on summary/smoothed ACF plot - figure(96), clf - subplot(2,1,1) - t=dt*(1:length(savedInputSignal)); - plot(t,savedInputSignal) - xlim([0 t(end)]) - title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']); - - - % compute ACF - switch saveAN_spikesOrProbability - case 'probability' - inputToACF=ANprobRateOutput.^0.5; - otherwise - inputToACF=ANoutput; - end - - disp ('computing ACF...') - [P, BFlist, sacf, boundaryValue] = ... - filteredSACF(inputToACF, method, filteredSACFParams); - disp(' ACF done.') - - % SACF - subplot(2,1,2) - imagesc(P) - ylabel('periodicities (Hz)') - xlabel('time (s)') - title(['running smoothed (root) SACF. ' saveAN_spikesOrProbability ' input']) - pt=[1 get(gca,'ytick')]; % force top xtick to show - set(gca,'ytick',pt) - set(gca,'ytickLabel', round(pitches(pt))) - tt=get(gca,'xtick'); - set(gca,'xtickLabel', round(100*t(tt))/100) -end - -path(restorePath) + + path(restorePath)
--- a/testPrograms/test_MAP1_14.m Tue Jun 07 10:54:42 2011 +0100 +++ b/testPrograms/test_MAP1_14.m Tue Jun 07 17:50:05 2011 +0100 @@ -9,7 +9,7 @@ % % #1 % Identify the file (in 'MAPparamsName') containing the model parameters -% +% % #2 % Identify the kind of model required (in 'AN_spikesOrProbability'). % A full brainstem model (spikes) can be computed or a shorter model @@ -49,7 +49,7 @@ % duration=0.020; % seconds sampleRate= 64000; % toneFrequency= 250:250:8000; % harmonic sequence (Hz) -toneFrequency= 2000; % or a pure tone (Hz8 +toneFrequency= 1000; % or a pure tone (Hz8 rampDuration=.005; % seconds @@ -58,25 +58,25 @@ fileName='twister_44kHz'; % fileName='new-da-44khz'; -% mix with an optional second file? +% ? and mix with an optional second file? mixerFile=[]; %or mixerFile='babble'; -leveldBSPL2=-60; +leveldBSPL2=60; %% #4 rms level % signal details -leveldBSPL= 60; % dB SPL +leveldBSPL= 70; % dB SPL %% #5 number of channels in the model % 21-channel model (log spacing) numChannels=21; -lowestBF=250; highestBF= 8000; +lowestBF=250; highestBF= 8000; BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels)); % or specify your own channel BFs -BFlist=toneFrequency; +% BFlist=toneFrequency; %% #6 change model parameters @@ -89,8 +89,7 @@ % constant of 80e-6 s (HSR fiber) when the probability option is selected. % paramChanges={'AN_IHCsynapseParams.ANspeedUpFactor=5;', ... % 'IHCpreSynapseParams.tauCa=86e-6;'}; -% paramChanges={'AN_IHCsynapseParams.ANspeedUpFactor=5;', ... -% 'DRNLParams.rateToAttenuationFactorProb = 0;'}; +% paramChanges={'DRNLParams.rateToAttenuationFactorProb = 0;'}; %% delare showMap options showMapOptions=[]; % use defaults @@ -120,14 +119,22 @@ case 'file' %% file input simple or mixed - [inputSignal sampleRate]=wavread(fileName); + [inputSignal sampleRate]=wavread(fileName); + dt=1/sampleRate; inputSignal=inputSignal(:,1); targetRMS=20e-6*10^(leveldBSPL/20); rms=(mean(inputSignal.^2))^0.5; amp=targetRMS/rms; inputSignal=inputSignal*amp; + silence= zeros(1,round(0.1/dt)); + inputSignal= [silence inputSignal' silence]; + if ~isempty(mixerFile) - [inputSignal2 sampleRate]=wavread(mixerFile); + [inputSignal2 sampleRate2]=wavread(mixerFile); + if ~isequal(sampleRate,sampleRate2) + error... + ('file and mixer file have different sample rates') + end inputSignal2=inputSignal2(:,1); [r c]=size(inputSignal); inputSignal2=inputSignal2(1:r); @@ -135,8 +142,18 @@ rms=(mean(inputSignal2.^2))^0.5; amp=targetRMS/rms; inputSignal2=inputSignal2*amp; + rampDuration=dt*length(inputSignal2)/3; + if rampDuration>0.5*duration, rampDuration=duration/2; end + rampTime=dt:dt:rampDuration; + time=dt*(1:length(inputSignal2)); + ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ... + ones(1,length(time)-length(rampTime))]; + inputSignal2=inputSignal2.*ramp'; + ramp=fliplr(ramp); + inputSignal2=inputSignal2.*ramp'; + end + inputSignal=inputSignal+inputSignal2; - end end @@ -159,7 +176,7 @@ % the model run is now complete. Now display the results showMAP(showMapOptions) for i=1:length(paramChanges) -disp(paramChanges{i}) + disp(paramChanges{i}) end toc
--- a/utilities/UTIL_makePSTH.m Tue Jun 07 10:54:42 2011 +0100 +++ b/utilities/UTIL_makePSTH.m Tue Jun 07 17:50:05 2011 +0100 @@ -1,4 +1,5 @@ -function PSTH =UTIL_makePSTH(inputData, dt, PSTHbinWidth) +function [PSTH dataPointsPerBin]= ... + UTIL_makePSTH(inputData, dt, PSTHbinWidth) % UTIL_PSTHmaker accumulates values into bins. % No corrections are applied % usage: