Revision 38:c2204b18f4a2 utilities
| utilities/UTIL_Butterworth.m | ||
|---|---|---|
| 1 |
function x=UTIL_Butterworth (x, dt, fl, fu, order) |
|
| 2 |
% UTIL_Butterworth (x, dt, fu, fl, order) |
|
| 3 |
% Taken from Yuel and beauchamp page 261 |
|
| 4 |
% NB error in their table for K (see their text) |
|
| 5 |
% x is original signal |
|
| 6 |
% x is the filtered output (approx 3 dB down at cutoff for first order filter) |
|
| 7 |
% fu, fl upper and lower cutoff |
|
| 8 |
% order is the number of times the filter is applied |
|
| 9 |
% (approx 6 dB attenuation is corrected) |
|
| 10 |
|
|
| 11 |
sampleRate=1/dt; |
|
| 12 |
if 4*fu>sampleRate |
|
| 13 |
error(['UTIL_Butterworth: sample rate ' num2str(sampleRate) ' is too low. Sampling rate should be >' num2str(4*fu)]) |
|
| 14 |
end |
|
| 15 |
|
|
| 16 |
|
|
| 17 |
q=(pi*dt*(fu-fl)); |
|
| 18 |
J=1/(1+ cot(q)); |
|
| 19 |
K= (2*cos(pi*dt*(fu+fl)))/((1+tan(q))*cos(q)); |
|
| 20 |
L= (tan(q)-1)/(tan(q)+1); |
|
| 21 |
b=[J 0 -J]; |
|
| 22 |
a=[1 -K -L]; |
|
| 23 |
|
|
| 24 |
[numChannels nDataPoints]=size(x); |
|
| 25 |
for channelNumber=1:numChannels |
|
| 26 |
for j=1:order |
|
| 27 |
% x2 compensate for 6 dB attenuation |
|
| 28 |
x(channelNumber,:)=2*filter(b, a, x(channelNumber,:)); |
|
| 29 |
end |
|
| 30 |
end |
|
| 31 |
|
|
| utilities/UTIL_CAPgenerator.m | ||
|---|---|---|
| 1 |
function wholeNerveCAP = UTIL_CAPgenerator... |
|
| 2 |
(ANresponse, dt, BFlist, numANfibers, plotCAP) |
|
| 3 |
% |
|
| 4 |
% Generates a compound action potential by convolving an impulse repsonse, |
|
| 5 |
% as defined by Mark Chertoff (2004, JASA), with the response of the |
|
| 6 |
% auditory nerve. |
|
| 7 |
% |
|
| 8 |
% |
|
| 9 |
% -e(-k*time)*SIN(2*PI()*f*time) |
|
| 10 |
% |
|
| 11 |
% mu(t) = e^-kt * sin(omega*t) |
|
| 12 |
% omega = 2 * pi * f |
|
| 13 |
% |
|
| 14 |
% |
|
| 15 |
% Robert T. Ferry |
|
| 16 |
% 01st May 2008 |
|
| 17 |
% |
|
| 18 |
% |
|
| 19 |
|
|
| 20 |
nChannels=length(BFlist); |
|
| 21 |
[r nSpikeEpochs]=size(ANresponse); |
|
| 22 |
|
|
| 23 |
wholeNerveCAP = []; |
|
| 24 |
channelCAP = []; |
|
| 25 |
e = exp(1); |
|
| 26 |
k = 1000; |
|
| 27 |
impulseDuration = 0.01; |
|
| 28 |
impulseFrequency = 750; |
|
| 29 |
impulseTime = dt:dt:impulseDuration; |
|
| 30 |
impulseResponse = -e.^(-k*impulseTime).*sin(2*pi*impulseFrequency*impulseTime); |
|
| 31 |
impulseResponse=impulseResponse-mean(impulseResponse); |
|
| 32 |
|
|
| 33 |
% WholeNerveCAP |
|
| 34 |
ANoutput = sum(ANresponse, 1); |
|
| 35 |
convolvedWholeNerveCAP = conv(ANoutput, impulseResponse(1,:)); |
|
| 36 |
% truncate |
|
| 37 |
convolvedWholeNerveCAP=convolvedWholeNerveCAP(1:nSpikeEpochs); |
|
| 38 |
|
|
| 39 |
% apply measurement time constant |
|
| 40 |
sampleRate=1/dt; |
|
| 41 |
upperFreq=sampleRate/4; |
|
| 42 |
lowPassCutOff=40; |
|
| 43 |
wholeNerveCAP = UTIL_Butterworth(convolvedWholeNerveCAP, dt, lowPassCutOff, upperFreq, 1); |
|
| 44 |
% or do not do this |
|
| 45 |
% wholeNerveCAP = convolvedWholeNerveCAP; |
|
| 46 |
|
|
| 47 |
% Plot output? |
|
| 48 |
|
|
| 49 |
CAPtime = dt:dt:dt*length(wholeNerveCAP); |
|
| 50 |
|
|
| 51 |
if (plotCAP == 1) |
|
| 52 |
figure(9) |
|
| 53 |
|
|
| 54 |
subplot(3,1,1), plot(impulseTime, impulseResponse) |
|
| 55 |
title('Impulse response')
|
|
| 56 |
xlabel('Time (s)'), ylabel('Amplitude (Pa)')
|
|
| 57 |
xlim([0 max(impulseTime)]), ylim([-inf inf]) |
|
| 58 |
|
|
| 59 |
subplot(3,1,2), plot(CAPtime, wholeNerveCAP) |
|
| 60 |
title(['AN CAP (whole-nerve) ' num2str(length(BFlist)) ' channels' num2str(numANfibers) ' fibers/ch']) |
|
| 61 |
xlabel('Time (s)'), ylabel('Amplitude (Pa)')
|
|
| 62 |
xlim([0 max(CAPtime)]), ylim([-inf inf]) |
|
| 63 |
|
|
| 64 |
end |
|
| utilities/UTIL_FFT.m | ||
|---|---|---|
| 1 |
function [fft_powerdB, fft_phase, frequencies, fft_ampdB]=... |
|
| 2 |
UTIL_FFT(getFFTsignal, dt, referenceAmplitude) |
|
| 3 |
% UTIL_FFT |
|
| 4 |
% [fft_powerdB, fft_phase, frequencies, fft_ampdB]= UTIL_FFT(getFFTsignal, dt, referenceAmplitude) |
|
| 5 |
|
|
| 6 |
x=length(getFFTsignal); % adjust the length of the signal to a power of 2 for FFT |
|
| 7 |
n=2; |
|
| 8 |
while n<=x |
|
| 9 |
n=n*2; |
|
| 10 |
end |
|
| 11 |
n=round(n/2); |
|
| 12 |
getFFTsignal=getFFTsignal(1:n); |
|
| 13 |
frequencies = (1/dt)*(1:n/2)/n; |
|
| 14 |
|
|
| 15 |
fft_result = fft(getFFTsignal, n); % Compute FFT of the input signal. |
|
| 16 |
fft_power = fft_result .* conj(fft_result); % / n; % Compute power spectrum. |
|
| 17 |
% Dividing by 'n' we get the power spectral density. |
|
| 18 |
% fft_power = fft_result .* conj(fft_result) / n; % Compute power spectralDensity. |
|
| 19 |
|
|
| 20 |
% Prepare the FFT for display |
|
| 21 |
fft_power=fft_power(1:length(fft_power)/2); % remove mirror frequencies |
|
| 22 |
% fft_powerdB = UTIL_amp2dB (fft_power, max(fft_power)); % convert to dB |
|
| 23 |
fft_powerdB = 10*log10(fft_power); % convert to dB |
|
| 24 |
|
|
| 25 |
% amplitude spectrum |
|
| 26 |
if nargin<3 |
|
| 27 |
referenceAmplitude=28e-6; % for SPL |
|
| 28 |
end |
|
| 29 |
% refAmpdB= referenceAmplitude^0.5; |
|
| 30 |
fft_ampdB = 10*log10(fft_power/referenceAmplitude); % convert to dB |
|
| 31 |
% peak_fft_powerdBSPL=max(fft_ampdB) |
|
| 32 |
|
|
| 33 |
|
|
| 34 |
fft_phase = angle(fft_result); % Compute the phase spectrum. |
|
| 35 |
fft_phase=fft_phase(1:length(fft_phase)/2); % remove mirror phases |
|
| 36 |
jags=find(diff(fft_phase)>0); % unwrap phase |
|
| 37 |
for i=jags, |
|
| 38 |
fft_phase(i+1:end)=fft_phase(i+1:end)-2*pi; |
|
| 39 |
end |
|
| utilities/UTIL_cascadePlot.m | ||
|---|---|---|
| 1 |
function UTIL_cascadePlot(toPlot, colValues) |
|
| 2 |
% % useful code |
|
| 3 |
[nChannels nLags]=size(toPlot); |
|
| 4 |
|
|
| 5 |
% cunning code to represent channels as parallel lines |
|
| 6 |
[nRows nCols]=size(toPlot); |
|
| 7 |
if nChannels<2 |
|
| 8 |
error('UTIL_cascadePlot: only one row found')
|
|
| 9 |
end |
|
| 10 |
|
|
| 11 |
% a is the height to be added to each channel |
|
| 12 |
a=max(max(toPlot))*(0:nRows-1)'; |
|
| 13 |
|
|
| 14 |
% peakGain emphasises the peak height |
|
| 15 |
% peaks can be higher than the space between channels |
|
| 16 |
peakGain=10; |
|
| 17 |
x=peakGain*toPlot+repmat(a,1,nCols); |
|
| 18 |
x=nRows*x/max(max(x)); |
|
| 19 |
|
|
| 20 |
for row=1:nRows-1 |
|
| 21 |
x(row,:)=min(x(row,:), x(row+1,:)); |
|
| 22 |
end |
|
| 23 |
|
|
| 24 |
plot(colValues, x','k') |
|
| 25 |
ylim([0 nRows]) |
|
| 26 |
|
|
| utilities/UTIL_plotMatrix.m | ||
|---|---|---|
| 2 | 2 |
% UTIL_plotMatrix general purpose plotting utility for plotting the results |
| 3 | 3 |
% of the MAP auditory model. |
| 4 | 4 |
% All plots are placed in subplots of a figure (default figure 1). |
| 5 |
% |
|
| 5 | 6 |
% Input arguments: |
| 6 | 7 |
% 'toPlot' is matrix (either numeric or logical) |
| 7 | 8 |
% 'method' is a structure containing plot instructions |
| 8 | 9 |
% |
| 9 |
% surface plots have log z-axis when all values are >1 |
|
| 10 |
% |
|
| 11 |
% when calling this function, always increment the subPlotNo in method |
|
| 12 |
% method.subPlotNo=method.subPlotNo+1; |
|
| 13 |
% method.subPlotNo=method.subPlotNo; |
|
| 14 |
% |
|
| 15 | 10 |
% mandatory parameters: |
| 16 |
% method.displaydt xValues spacing between data points |
|
| 17 |
% method.numPlots total number of subPlots in the figure |
|
| 18 |
% method.subPlotNo number of this plot |
|
| 19 |
% method.yValues mandatory only for 3D plots |
|
| 11 |
% method.displaydt xValues spacing between data points |
|
| 12 |
% method.yValues yaxis labels mandatory only for 3D plots |
|
| 20 | 13 |
% |
| 21 | 14 |
% optional |
| 22 |
% method.figureNo normally figure(1) |
|
| 23 |
% method.zValuesRange [min max] value pair to define yaxis limits |
|
| 24 |
% method.zValuesRange [min max] CLIMS for 3-D plot |
|
| 25 |
% method.yLabel (string) y-axis label |
|
| 26 |
% method.xLabel (string) x-axis label |
|
| 27 |
% method.title (string) subplot title |
|
| 15 |
% method.figureNo default figure(1) |
|
| 16 |
% method.numPlots number of subPlots in the figure (default=1) |
|
| 17 |
% method.subPlotNo number of this plot (default=1) |
|
| 18 |
% method.zValuesRange [min max] value pair to define yaxis limits |
|
| 19 |
% method.zValuesRange [min max] CLIMS for 3-D plot |
|
| 20 |
% method.yLabel (string) y-axis label |
|
| 21 |
% method.minyMaxy y-axis limits |
|
| 22 |
% method.xLabel (string) x-axis label |
|
| 23 |
% method.title (string) subplot title |
|
| 28 | 24 |
% method.bar =1, to force bar histogram (single channel only) |
| 29 |
% method.view 3D plot 'view' settings e.g. [-6 40] |
|
| 30 |
% method.axes (handle) used for writing to GUIs (specifies panel) |
|
| 31 |
% method.maxPixels maximum number of pixels (used to speed plotting) |
|
| 32 |
% method.blackOnWhite =1; % NB inverts display for 2D plots |
|
| 33 |
% method.forceLog positive values are put on log z-scale |
|
| 34 |
% method.rasterDotSize min value is 1 |
|
| 35 |
% method.defaultFontSize |
|
| 36 |
% method.timeStart default=dt |
|
| 37 |
% method.defaultTextColor default ='w' |
|
| 38 |
% method.defaultAxesColor default = 'w' |
|
| 39 |
% method.nCols default=1 |
|
| 40 |
% method.nRows default=method.numPlots |
|
| 25 |
% method.view 3D plot 'view' settings e.g. [-6 40] |
|
| 26 |
% method.axes (handle) where to plot (overules all others) |
|
| 27 |
% method.maxPixels maximum number of pixels (used to speed plotting) |
|
| 28 |
% method.blackOnWhite =1; inverts display for 2D plots |
|
| 29 |
% method.forceLog positive values are put on log z-scale |
|
| 30 |
% method.rasterDotSize min value is 1 |
|
| 31 |
% method.defaultFontSize deafult= 12 |
|
| 32 |
% method.timeStart default= dt |
|
| 33 |
% method.defaultTextColor default ='k' |
|
| 34 |
% method.defaultAxesColor default ='k' |
|
| 35 |
% method.nCols default = 1 (layout for subplots) |
|
| 36 |
% method.nRows default=method.numPlots (layout for subplots |
|
| 37 |
% method.segmentNumber plot only this segment while 'hold on' |
|
| 41 | 38 |
% |
| 42 |
% useful paste for calling program |
|
| 43 |
% method.numPlots=method.numPlots; |
|
| 44 |
% method.subPlotNo=method.subPlotNo+1; |
|
| 45 |
% method.subPlotNo=method.subPlotNo; |
|
| 46 |
% dt=dt; |
|
| 47 |
% method.yValues=method.nonlinCF; % for 3D plots only |
|
| 48 |
% |
|
| 49 |
% method.figureNo=1; |
|
| 50 |
% method.yLabel='useThisLabel'; |
|
| 51 |
% method.xLabel='use this label'; |
|
| 52 |
% method.title='myTitle'; |
|
| 53 |
% |
|
| 39 |
% e.g. |
|
| 54 | 40 |
% UTIL_plotMatrix(toPlot, method) |
| 55 | 41 |
|
| 56 | 42 |
dt=method.displaydt; |
| 43 |
[r cols]=size(toPlot); |
|
| 44 |
if cols==1 |
|
| 45 |
% toPlot should be a wide matrix or a long vector |
|
| 46 |
toPlot=toPlot'; |
|
| 47 |
end |
|
| 48 |
|
|
| 49 |
if ~isfield(method,'numPlots') || isempty(method.numPlots) |
|
| 50 |
method.numPlots =1; |
|
| 51 |
method.subPlotNo =1; |
|
| 52 |
end |
|
| 53 |
|
|
| 57 | 54 |
if ~isfield(method,'figureNo') || isempty(method.figureNo) |
| 58 | 55 |
method.figureNo=99; |
| 59 | 56 |
end |
| 57 |
|
|
| 60 | 58 |
% if ~isfield(method,'zValuesRange') || isempty(method.zValuesRange) |
| 61 | 59 |
% method.zValuesRange=[-inf inf]; |
| 62 | 60 |
% end |
| 63 | 61 |
|
| 64 | 62 |
% set some defaults |
| 65 |
if ~isfield( method,'plotDivider') || isempty(method.plotDivider) |
|
| 66 |
method.plotDivider=0; |
|
| 67 |
end |
|
| 68 | 63 |
if ~isfield( method,'blackOnWhite') || isempty(method.blackOnWhite) |
| 69 |
method.blackOnWhite=0;
|
|
| 64 |
method.blackOnWhite=0; |
|
| 70 | 65 |
end |
| 71 | 66 |
if ~isfield(method,'timeStart')|| isempty(method.timeStart) |
| 72 |
method.timeStart=dt;
|
|
| 67 |
method.timeStart=dt; |
|
| 73 | 68 |
end |
| 74 | 69 |
if ~isfield(method,'objectDuration') || isempty(method.objectDuration) |
| 75 | 70 |
[nRows nCols]=size(toPlot); method.objectDuration=dt*nCols; |
| 76 | 71 |
end |
| 77 | 72 |
if ~isfield(method,'defaultFontSize') || isempty(method.defaultFontSize) |
| 78 |
method.defaultFontSize=12;
|
|
| 73 |
method.defaultFontSize=12; |
|
| 79 | 74 |
end |
| 80 | 75 |
if ~isfield(method,'defaultTextColor') || isempty(method.defaultTextColor) |
| 81 | 76 |
method.defaultTextColor='k'; |
| ... | ... | |
| 84 | 79 |
defaultTextColor='k'; |
| 85 | 80 |
end |
| 86 | 81 |
if ~isfield( method,'defaultAxesColor') || isempty(method.defaultAxesColor) |
| 87 |
method.defaultAxesColor=defaultTextColor;
|
|
| 82 |
method.defaultAxesColor=defaultTextColor; |
|
| 88 | 83 |
end |
| 89 | 84 |
defaultAxesColor=method.defaultAxesColor; |
| 90 | 85 |
|
| ... | ... | |
| 107 | 102 |
% if both are specified, 'axes' takes priority |
| 108 | 103 |
figure(method.figureNo) |
| 109 | 104 |
if isfield(method,'axes') && ~isempty(method.axes) |
| 110 |
% select an axis in some location other than 'figure' |
|
| 111 |
h=axes(method.axes); |
|
| 105 |
% user defines where to plot it |
|
| 106 |
axes(method.axes); |
|
| 107 |
method.numPlots =1; |
|
| 108 |
method.subPlotNo =1; |
|
| 109 |
|
|
| 112 | 110 |
else |
| 113 | 111 |
% now using a regular figure |
| 114 | 112 |
if method.subPlotNo>method.numPlots; |
| ... | ... | |
| 119 | 117 |
|
| 120 | 118 |
if isfield(method,'segmentNumber') && ~isempty(method.segmentNumber)... |
| 121 | 119 |
&& method.segmentNumber>1 |
| 122 |
% in multi-segment mode do not clear the image
|
|
| 120 |
% in multi-segment mode do not clear the image |
|
| 123 | 121 |
% from the previous segment |
| 124 | 122 |
hold on |
| 125 | 123 |
else |
| ... | ... | |
| 139 | 137 |
yValues=1:numYvalues; |
| 140 | 138 |
end |
| 141 | 139 |
|
| 140 |
if round(numYvalues/length(yValues))>1 |
|
| 141 |
% case where the plot matrix is double height (e.g. LSR+HSR) |
|
| 142 |
yValues=[yValues yValues]; |
|
| 143 |
method.plotDivider=1; |
|
| 144 |
else |
|
| 145 |
method.plotDivider=0; |
|
| 146 |
end |
|
| 147 |
|
|
| 142 | 148 |
% Now start the plot. |
| 143 | 149 |
% 3D plotting for 4 or more channels |
| 144 | 150 |
% otherwise special cases for fewer channels |
| ... | ... | |
| 195 | 201 |
end |
| 196 | 202 |
|
| 197 | 203 |
otherwise % >3 channels: surface plot |
| 198 |
% add white line to separate HSR and LSR
|
|
| 204 |
% add line to separate HSR and LSR |
|
| 199 | 205 |
if method.plotDivider && size(toPlot,1) > 2 |
| 200 | 206 |
[r c]=size(toPlot); |
| 201 |
% if isempty(method.zValuesRange), method.zValuesRange=0; end |
|
| 202 |
% mm=method.zValuesRange(2); |
|
| 203 | 207 |
emptyLine=max(max(toPlot))*ones(2,c); |
| 204 | 208 |
halfway=round(r/2); |
| 205 | 209 |
toPlot=[toPlot(1:halfway,:); emptyLine; toPlot(halfway+1:end,:)]; |
| ... | ... | |
| 225 | 229 |
% zValuesRange |
| 226 | 230 |
if isfield(method,'zValuesRange') ... |
| 227 | 231 |
&& ~isempty(method.zValuesRange) |
| 228 |
% zValuesRange gives the max and min values |
|
| 229 |
% a=method.zValuesRange(1); |
|
| 230 |
% b=method.zValuesRange(2); |
|
| 231 |
% toPlot=(toPlot-a)/(b-a); |
|
| 232 |
% clims=[0 1]; |
|
| 233 | 232 |
clims=(method.zValuesRange); |
| 234 | 233 |
imagesc(xValues, yValues, toPlot, clims), axis xy; %NB assumes equally spaced y-values |
| 235 | 234 |
else |
| 236 | 235 |
% automatically scaled |
| 237 | 236 |
imagesc(xValues, yValues, toPlot), axis xy; %NB assumes equally spaced y-values |
| 238 | 237 |
|
| 239 |
% if ~isfield(method,'zValuesRange') |
|
| 240 |
% method.zValuesRange=[-inf inf]; |
|
| 241 |
% end |
|
| 242 |
% |
|
| 243 |
% if method.blackOnWhite |
|
| 244 |
% % NB plotted values have negative sign for black on white |
|
| 245 |
% caxis([-method.zValuesRange(2) -method.zValuesRange(1)]) |
|
| 246 |
% else |
|
| 247 |
% caxis(method.zValuesRange) |
|
| 248 |
% end |
|
| 249 |
|
|
| 250 | 238 |
if ~isfield(method,'zValuesRange')... |
| 251 | 239 |
|| isempty(method.zValuesRange) |
| 252 | 240 |
method.zValuesRange=[-inf inf]; |
| ... | ... | |
| 264 | 252 |
% NB segmentation may shorten signal duration |
| 265 | 253 |
[r c]=size(toPlot); |
| 266 | 254 |
imageDuration=c*method.displaydt; |
| 267 |
xlim([0 imageDuration]) |
|
| 268 |
% xlim([0 method.objectDuration]) |
|
| 269 |
|
|
| 255 |
xlim([0 imageDuration]) |
|
| 256 |
|
|
| 270 | 257 |
% yaxis |
| 271 | 258 |
if isfield(method,'minyMaxy') && ~isempty(method.minyMaxy) |
| 272 | 259 |
ylim(method.minyMaxy) |
| ... | ... | |
| 275 | 262 |
ylim([min(yValues) max(yValues)]) |
| 276 | 263 |
end |
| 277 | 264 |
end |
| 278 |
if min(yValues)>1 % put channel array on a log scale |
|
| 265 |
|
|
| 266 |
% y-axis design yTickLabels |
|
| 267 |
if min(yValues)>1 |
|
| 279 | 268 |
tickValues=[min(yValues) max(yValues)]; |
| 269 |
tickLabels=num2str(tickValues'); |
|
| 270 |
if method.plotDivider && size(toPlot,1) > 2 |
|
| 271 |
% show min/max yvalues with slight shift |
|
| 272 |
yList=yValues; |
|
| 273 |
yValues=1:length(yValues); |
|
| 274 |
tickValues=[1 halfway-1 halfway+2 length(yValues)]; |
|
| 275 |
idx=[1 halfway halfway+1 length(yValues)]; |
|
| 276 |
tickLabels=num2str(yList(idx)'); |
|
| 277 |
imagesc(xValues, yValues, toPlot), axis xy; |
|
| 278 |
end |
|
| 279 |
|
|
| 280 | 280 |
set(gca,'ytick',tickValues) |
| 281 |
set(gca,'ytickLabel', strvcat(num2str(tickValues')))
|
|
| 281 |
set(gca,'ytickLabel', strvcat(tickLabels))
|
|
| 282 | 282 |
set(gca,'FontSize', method.defaultFontSize) |
| 283 | 283 |
end |
| 284 | 284 |
|
| 285 | 285 |
end |
| 286 | 286 |
|
| 287 |
else % is logical |
|
| 288 |
|
|
| 287 |
else % is logical |
|
| 289 | 288 |
% logical implies spike array. Use raster plot |
| 290 | 289 |
[y,x]=find(toPlot); %locate all spikes: y is fiber number ie row |
| 291 | 290 |
x=x*dt+method.timeStart; % x is time |
| ... | ... | |
| 314 | 313 |
plot([0 c*method.displaydt],[halfWayUp halfWayUp], 'b') |
| 315 | 314 |
hold off |
| 316 | 315 |
end |
| 316 |
|
|
| 317 | 317 |
end |
| 318 | 318 |
|
| 319 | 319 |
set(gca, 'xcolor', defaultAxesColor) |
| ... | ... | |
| 330 | 330 |
set(gca,'xtick',[],'FontSize', method.defaultFontSize) |
| 331 | 331 |
if method.subPlotNo==method.numPlots |
| 332 | 332 |
if isfield(method,'xLabel') && ~isempty(method.xLabel) |
| 333 |
% set(gca,'ActivePositionProperty','outerposition') |
|
| 333 |
% set(gca,'ActivePositionProperty','outerposition')
|
|
| 334 | 334 |
% xlabel(method.xLabel) |
| 335 | 335 |
xlabel(method.xLabel, 'FontSize', method.defaultFontSize, 'color', defaultTextColor) |
| 336 | 336 |
end |
| utilities/UTIL_printTabTable.m | ||
|---|---|---|
| 2 | 2 |
% printTabTable prints a matrix as a table with tabs |
| 3 | 3 |
%headers are optional |
| 4 | 4 |
%headers=strvcat('firstname', 'secondname')
|
| 5 |
% printTabTable([1 2; 3 4],strvcat('a1','a2'));
|
|
| 5 |
% UTIL_printTabTable([1 2; 3 4],strvcat('a1','a2'));
|
|
| 6 | 6 |
|
| 7 | 7 |
if nargin<3 |
| 8 | 8 |
format='%g'; |
| utilities/UTIL_showMAP.m | ||
|---|---|---|
| 1 |
function UTIL_showMAP (showMapOptions, paramChanges)
|
|
| 1 |
function UTIL_showMAP (showMapOptions) |
|
| 2 | 2 |
% UTIL_showMAP produces summaries of the output from MAP's mmost recent run |
| 3 |
% All MAP outputs are stored in global variables and UTIL_showMAP |
|
| 4 |
% simply assumes that they are in place.
|
|
| 3 |
% All MAP_14 outputs are stored in global variables and UTIL_showMAP
|
|
| 4 |
% simply assumes that they are in place. It uses whatever it there.
|
|
| 5 | 5 |
% |
| 6 |
% showMapOptions |
|
| 6 |
% showMapOptions:
|
|
| 7 | 7 |
% showMapOptions.printModelParameters=1; % print model parameters |
| 8 | 8 |
% showMapOptions.showModelOutput=1; % plot all stages output |
| 9 | 9 |
% showMapOptions.printFiringRates=1; % mean activity at all stages |
| 10 | 10 |
% showMapOptions.showACF=1; % SACF (probabilities only) |
| 11 | 11 |
% showMapOptions.showEfferent=1; % plot of efferent activity |
| 12 |
% showMapOptions.surfProbability=0; % HSR (probability) surf plot |
|
| 13 |
% showMapOptions.fileName=[]; % parameter filename for plot title |
|
| 14 |
|
|
| 12 |
% showMapOptions.surfAN=0; % AN output surf plot |
|
| 13 |
% showMapOptions.fileName=''; % parameter filename for plot title |
|
| 14 |
% showMapOptions.PSTHbinwidth=0.001 % binwidth for surface plots |
|
| 15 |
% showMapOptions.colorbar=1; % request color bar if appropriate |
|
| 16 |
% showMapOptions.view=[0 90]; |
|
| 15 | 17 |
% dbstop if warning |
| 16 | 18 |
|
| 17 |
global dt ANdt savedBFlist saveAN_spikesOrProbability saveMAPparamsName... |
|
| 18 |
savedInputSignal OMEextEarPressure TMoutput OMEoutput ARattenuation ... |
|
| 19 |
DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV... |
|
| 20 |
IHCoutput ANprobRateOutput ANoutput savePavailable ANtauCas ... |
|
| 21 |
CNtauGk CNoutput ICoutput ICmembraneOutput ICfiberTypeRates ... |
|
| 22 |
MOCattenuation |
|
| 23 |
global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams |
|
| 24 |
global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams |
|
| 19 |
% Discover results left behind by MAP1_14 |
|
| 20 |
global savedParamChanges savedBFlist saveAN_spikesOrProbability ... |
|
| 21 |
saveMAPparamsName savedInputSignal dt dtSpikes ... |
|
| 22 |
OMEextEarPressure TMoutput OMEoutput DRNLoutput... |
|
| 23 |
IHC_cilia_output IHCrestingCiliaCond IHCrestingV... |
|
| 24 |
IHCoutput ANprobRateOutput ANoutput savePavailable saveNavailable ... |
|
| 25 |
ANtauCas CNtauGk CNoutput ICoutput ICmembraneOutput ICfiberTypeRates... |
|
| 26 |
MOCattenuation ARattenuation |
|
| 27 |
|
|
| 28 |
% Find parameters created in MAPparams<name> |
|
| 29 |
global inputStimulusParams OMEParams DRNLParams IHC_cilia_RPParams |
|
| 30 |
global IHCpreSynapseParams AN_IHCsynapseParams |
|
| 31 |
global MacGregorParams MacGregorMultiParams filteredSACFParams |
|
| 25 | 32 |
global ICrate |
| 26 | 33 |
|
| 27 | 34 |
|
| 28 | 35 |
restorePath=path; |
| 29 | 36 |
addpath ( ['..' filesep 'utilities'], ['..' filesep 'parameterStore']) |
| 30 | 37 |
|
| 31 |
if nargin<1 |
|
| 32 |
showMapOptions=[];
|
|
| 33 |
end |
|
| 38 |
if nargin<1, showMapOptions=[]; end
|
|
| 39 |
paramChanges=savedParamChanges;
|
|
| 40 |
|
|
| 34 | 41 |
% defaults (plot staged outputs and print rates only) if options omitted |
| 35 | 42 |
if ~isfield(showMapOptions,'printModelParameters') |
| 36 | 43 |
showMapOptions.printModelParameters=0; end |
| 37 |
if ~isfield(showMapOptions,'showModelOutput'),showMapOptions.showModelOutput=1;end |
|
| 38 |
if ~isfield(showMapOptions,'printFiringRates'),showMapOptions.printFiringRates=1;end |
|
| 44 |
if ~isfield(showMapOptions,'showModelOutput'),showMapOptions.showModelOutput=0;end |
|
| 45 |
if ~isfield(showMapOptions,'printFiringRates'),showMapOptions.printFiringRates=0;end |
|
| 46 |
if ~isfield(showMapOptions,'surfAN'),showMapOptions.surfAN=0;end |
|
| 47 |
if ~isfield(showMapOptions,'ICrates'),showMapOptions.ICrates=0;end |
|
| 39 | 48 |
if ~isfield(showMapOptions,'showACF'),showMapOptions.showACF=0;end |
| 40 | 49 |
if ~isfield(showMapOptions,'showEfferent'),showMapOptions.showEfferent=0;end |
| 41 |
if ~isfield(showMapOptions,'surfProbability'),showMapOptions.surfProbability=0;end
|
|
| 50 |
if ~isfield(showMapOptions,'PSTHbinwidth'),showMapOptions.PSTHbinwidth=0.001;end
|
|
| 42 | 51 |
if ~isfield(showMapOptions,'fileName'),showMapOptions.fileName=[];end |
| 43 |
if ~isfield(showMapOptions,'surfSpikes'),showMapOptions.surfSpikes=0;end |
|
| 44 |
if ~isfield(showMapOptions,'ICrates'),showMapOptions.ICrates=0;end |
|
| 52 |
if ~isfield(showMapOptions,'colorbar'),showMapOptions.colorbar=1;end |
|
| 53 |
if ~isfield(showMapOptions,'view'),showMapOptions.view=[-25 80];end |
|
| 54 |
if ~isfield(showMapOptions,'ICrasterPlot'),showMapOptions.ICrasterPlot=0;end |
|
| 45 | 55 |
|
| 46 | 56 |
%% send all model parameters to command window |
| 47 | 57 |
if showMapOptions.printModelParameters |
| 48 | 58 |
% Read parameters from MAPparams<***> file in 'parameterStore' folder |
| 49 | 59 |
% and print out all parameters |
| 50 |
if nargin>1 |
|
| 51 | 60 |
cmd=['MAPparams' saveMAPparamsName '(-1, 1/dt, 1, paramChanges);']; |
| 52 | 61 |
eval(cmd); |
| 53 |
else |
|
| 54 |
cmd=['MAPparams' saveMAPparamsName '(-1, 1/dt, 1);']; |
|
| 55 |
eval(cmd); |
|
| 56 |
disp(' no parameter changes found')
|
|
| 57 |
end |
|
| 58 | 62 |
end |
| 59 | 63 |
|
| 64 |
% ignore zero elements in signal |
|
| 65 |
signalRMS=mean(savedInputSignal(find(savedInputSignal)).^2)^0.5; |
|
| 66 |
signalRMSdb=20*log10(signalRMS/20e-6); |
|
| 67 |
nANfiberTypes=length(ANtauCas); |
|
| 68 |
|
|
| 60 | 69 |
%% summarise firing rates in command window |
| 61 | 70 |
if showMapOptions.printFiringRates |
| 62 | 71 |
%% print summary firing rates |
| ... | ... | |
| 64 | 73 |
disp('summary')
|
| 65 | 74 |
disp(['AR: ' num2str(min(ARattenuation))]) |
| 66 | 75 |
disp(['MOC: ' num2str(min(min(MOCattenuation)))]) |
| 67 |
nANfiberTypes=length(ANtauCas); |
|
| 68 | 76 |
if strcmp(saveAN_spikesOrProbability, 'spikes') |
| 69 | 77 |
nANfibers=size(ANoutput,1); |
| 70 | 78 |
nHSRfibers=nANfibers/nANfiberTypes; |
| ... | ... | |
| 93 | 101 |
%% figure (99): display output from all model stages |
| 94 | 102 |
if showMapOptions.showModelOutput |
| 95 | 103 |
plotInstructions.figureNo=99; |
| 96 |
|
|
| 97 |
% ignore zero elements in signal |
|
| 98 |
signalRMS=mean(savedInputSignal(find(savedInputSignal)).^2)^0.5; |
|
| 99 |
signalRMSdb=20*log10(signalRMS/20e-6); |
|
| 100 | 104 |
|
| 101 | 105 |
% plot signal (1) |
| 102 | 106 |
plotInstructions.displaydt=dt; |
| ... | ... | |
| 104 | 108 |
plotInstructions.subPlotNo=1; |
| 105 | 109 |
plotInstructions.title=... |
| 106 | 110 |
['stimulus (Pa). ' num2str(signalRMSdb, '%4.0f') ' dB SPL']; |
| 107 |
r=size(savedInputSignal,1); |
|
| 108 |
if r==1, savedInputSignal=savedInputSignal'; end |
|
| 109 |
UTIL_plotMatrix(savedInputSignal', plotInstructions); |
|
| 111 |
UTIL_plotMatrix(savedInputSignal, plotInstructions); |
|
| 110 | 112 |
|
| 111 | 113 |
% stapes (2) |
| 112 | 114 |
plotInstructions.subPlotNo=2; |
| ... | ... | |
| 119 | 121 |
[r c]=size(DRNLoutput); |
| 120 | 122 |
if r>1 |
| 121 | 123 |
plotInstructions.title= ['BM displacement']; |
| 122 |
UTIL_plotMatrix(abs(DRNLoutput), plotInstructions); |
|
| 124 |
UTIL_plotMatrix(abs(DRNLoutput), plotInstructions);
|
|
| 123 | 125 |
else |
| 124 | 126 |
plotInstructions.title= ['BM displacement. BF=' ... |
| 125 | 127 |
num2str(savedBFlist) ' Hz']; |
| 126 |
UTIL_plotMatrix(DRNLoutput, plotInstructions); |
|
| 128 |
UTIL_plotMatrix(DRNLoutput, plotInstructions);
|
|
| 127 | 129 |
end |
| 128 | 130 |
|
| 129 | 131 |
switch saveAN_spikesOrProbability |
| 130 | 132 |
case 'spikes' |
| 131 | 133 |
% AN (4) |
| 132 |
plotInstructions.displaydt=ANdt;
|
|
| 134 |
plotInstructions.displaydt=dtSpikes;
|
|
| 133 | 135 |
plotInstructions.title='AN'; |
| 134 | 136 |
plotInstructions.subPlotNo=4; |
| 135 | 137 |
plotInstructions.yLabel='BF'; |
| ... | ... | |
| 146 | 148 |
UTIL_plotMatrix(ANoutput, plotInstructions); |
| 147 | 149 |
|
| 148 | 150 |
% CN (5) |
| 149 |
plotInstructions.displaydt=ANdt;
|
|
| 151 |
plotInstructions.displaydt=dtSpikes;
|
|
| 150 | 152 |
plotInstructions.subPlotNo=5; |
| 151 | 153 |
plotInstructions.title='CN spikes'; |
| 152 | 154 |
if sum(sum(CNoutput))<100 |
| ... | ... | |
| 155 | 157 |
UTIL_plotMatrix(CNoutput, plotInstructions); |
| 156 | 158 |
|
| 157 | 159 |
% IC (6) |
| 158 |
plotInstructions.displaydt=ANdt;
|
|
| 160 |
plotInstructions.displaydt=dtSpikes;
|
|
| 159 | 161 |
plotInstructions.subPlotNo=6; |
| 160 | 162 |
plotInstructions.title='Brainstem 2nd level'; |
| 161 | 163 |
if size(ICoutput,1)>1 |
| ... | ... | |
| 172 | 174 |
end |
| 173 | 175 |
|
| 174 | 176 |
otherwise % AN rate based on probability of firing |
| 175 |
PSTHbinWidth=0.001;
|
|
| 177 |
PSTHbinWidth=0.0005;
|
|
| 176 | 178 |
PSTH= UTIL_PSTHmakerb(ANprobRateOutput, dt, PSTHbinWidth); |
| 177 |
% PSTH = makeANsmooth(PSTH, 1/PSTHbinWidth);
|
|
| 179 |
plotInstructions=[];
|
|
| 178 | 180 |
plotInstructions.displaydt=PSTHbinWidth; |
| 181 |
plotInstructions.plotDivider=0; |
|
| 179 | 182 |
plotInstructions.numPlots=2; |
| 180 | 183 |
plotInstructions.subPlotNo=2; |
| 181 | 184 |
plotInstructions.yLabel='BF'; |
| 182 |
plotInstructions.xLabel='time'; |
|
| 185 |
plotInstructions.yValues=savedBFlist; |
|
| 186 |
plotInstructions.xLabel='time (s)'; |
|
| 183 | 187 |
plotInstructions.zValuesRange= [0 300]; |
| 188 |
|
|
| 184 | 189 |
if nANfiberTypes>1, |
| 185 |
plotInstructions.yLabel='LSR HSR'; |
|
| 190 |
plotInstructions.yLabel='LSR HSR';
|
|
| 186 | 191 |
plotInstructions.plotDivider=1; |
| 187 | 192 |
end |
| 188 | 193 |
plotInstructions.title='AN - spike rate'; |
| 189 | 194 |
UTIL_plotMatrix(PSTH, plotInstructions); |
| 190 | 195 |
shading interp |
| 191 |
colorbar('southOutside')
|
|
| 196 |
if showMapOptions.colorbar, colorbar('southOutside'), end
|
|
| 192 | 197 |
end |
| 193 | 198 |
set(gcf,'name','MAP output') |
| 194 | 199 |
end |
| 195 | 200 |
|
| 196 |
if showMapOptions.surfProbability &&... |
|
| 201 |
|
|
| 202 |
%% surface plot of AN response |
|
| 203 |
% probability |
|
| 204 |
if showMapOptions.surfAN &&... |
|
| 197 | 205 |
strcmp(saveAN_spikesOrProbability,'probability') && ... |
| 198 | 206 |
length(savedBFlist)>2 |
| 199 |
%% surface plot of probability |
|
| 200 |
% select only HSR fibers |
|
| 201 |
figure(97), clf |
|
| 202 |
HSRprobOutput= ANprobRateOutput(end-length(savedBFlist)+1:end,:); |
|
| 203 |
PSTHbinWidth=0.001; |
|
| 204 |
PSTH=UTIL_PSTHmakerb(HSRprobOutput, ANdt, PSTHbinWidth); |
|
| 205 |
[nY nX]=size(PSTH); |
|
| 206 |
time=PSTHbinWidth*(1:nX); |
|
| 207 |
surf(time, savedBFlist, PSTH) |
|
| 208 |
shading interp |
|
| 209 |
set(gca, 'yScale','log') |
|
| 210 |
xlim([0 max(time)]) |
|
| 211 |
ylim([0 max(savedBFlist)]) |
|
| 212 |
zlim([0 1000]) |
|
| 213 |
xlabel('time (s)')
|
|
| 214 |
ylabel('best frequency (Hz)')
|
|
| 215 |
zlabel('spike rate')
|
|
| 216 |
view([-20 60]) |
|
| 217 |
% view([0 90]) |
|
| 218 |
disp(['max max AN: ' num2str(max(max(PSTH)))]) |
|
| 219 |
|
|
| 220 |
title (['firing probability of HSR fibers only. Level= ' ... |
|
| 221 |
num2str(signalRMSdb,'% 3.0f') ' dB']) |
|
| 222 |
colorbar('southOutside')
|
|
| 207 |
% select only HSR fibers |
|
| 208 |
figure(97), clf |
|
| 209 |
PSTHbinWidth=showMapOptions.PSTHbinwidth; |
|
| 210 |
PSTH= UTIL_PSTHmakerb(... |
|
| 211 |
ANprobRateOutput(end-length(savedBFlist)+1:end,:), ... |
|
| 212 |
dt, PSTHbinWidth); |
|
| 213 |
[nY nX]=size(PSTH); |
|
| 214 |
time=PSTHbinWidth*(1:nX); |
|
| 215 |
surf(time, savedBFlist, PSTH) |
|
| 216 |
caxis([0 500]) |
|
| 217 |
shading interp |
|
| 218 |
set(gca, 'yScale','log') |
|
| 219 |
xlim([0 max(time)]) |
|
| 220 |
ylim([0 max(savedBFlist)]) |
|
| 221 |
zlim([0 500]) |
|
| 222 |
myFontSize=10; |
|
| 223 |
xlabel('time (s)', 'fontsize',myFontSize)
|
|
| 224 |
ylabel('BF (Hz)', 'fontsize',myFontSize)
|
|
| 225 |
zlabel('spike rate')
|
|
| 226 |
view(showMapOptions.view) |
|
| 227 |
set(gca,'ytick',[500 1000 2000 8000],'fontSize',myFontSize) |
|
| 228 |
title (['AN response. Level= ' ... |
|
| 229 |
num2str(signalRMSdb,'% 3.0f') ' dB'... |
|
| 230 |
' binwidth= ' num2str(1000*PSTHbinWidth) ' s']) |
|
| 231 |
if showMapOptions.colorbar, colorbar('southOutside'), end
|
|
| 223 | 232 |
end |
| 224 | 233 |
|
| 225 |
%% surface plot of AN spikes |
|
| 226 |
if showMapOptions.surfSpikes ... |
|
| 227 |
&& strcmp(saveAN_spikesOrProbability, 'spikes') |
|
| 234 |
%% surfAN ('spikes')
|
|
| 235 |
if showMapOptions.surfAN ... |
|
| 236 |
&& strcmp(saveAN_spikesOrProbability, 'spikes')... |
|
| 237 |
&& length(savedBFlist)>2 |
|
| 228 | 238 |
figure(97), clf |
| 239 |
% combine fibers across channels |
|
| 240 |
nFibersPerChannel=AN_IHCsynapseParams.numFibers; |
|
| 241 |
[r nEpochs]=size(ANoutput); |
|
| 242 |
nChannels=round(r/nFibersPerChannel); |
|
| 243 |
x=reshape(ANoutput,nChannels,nFibersPerChannel,nEpochs); |
|
| 244 |
x=squeeze(sum(x,2)); |
|
| 229 | 245 |
% select only HSR fibers at the bottom of the matrix |
| 230 |
ANoutput= ANoutput(end-length(savedBFlist)+1:end,:); |
|
| 231 |
PSTHbinWidth=0.005; % 1 ms bins |
|
| 232 |
PSTH=UTIL_makePSTH(ANoutput, ANdt, PSTHbinWidth); |
|
| 246 |
HSRoutput= x(end-length(savedBFlist)+1:end,:); |
|
| 247 |
PSTH=HSRoutput; |
|
| 248 |
PSTHbinWidth=showMapOptions.PSTHbinwidth; |
|
| 249 |
PSTH=UTIL_makePSTH(HSRoutput, dtSpikes, PSTHbinWidth); |
|
| 233 | 250 |
[nY nX]=size(PSTH); |
| 234 | 251 |
time=PSTHbinWidth*(1:nX); |
| 235 | 252 |
surf(time, savedBFlist, PSTH) |
| ... | ... | |
| 241 | 258 |
xlabel('time (s)')
|
| 242 | 259 |
ylabel('best frequency (Hz)')
|
| 243 | 260 |
zlabel('spike rate')
|
| 244 |
view([-20 60]) |
|
| 245 |
% view([0 90]) |
|
| 261 |
view(showMapOptions.view) |
|
| 246 | 262 |
title ([showMapOptions.fileName ': ' num2str(signalRMSdb,'% 3.0f') ' dB']) |
| 247 | 263 |
set(97,'name', 'spikes surface plot') |
| 248 | 264 |
end |
| 249 | 265 |
|
| 266 |
%% IC raster plot |
|
| 267 |
if showMapOptions.ICrasterPlot &&... |
|
| 268 |
strcmp(saveAN_spikesOrProbability,'spikes') && ... |
|
| 269 |
length(savedBFlist)>2 |
|
| 270 |
figure(91), clf |
|
| 271 |
plotInstructions=[]; |
|
| 272 |
plotInstructions.numPlots=2; |
|
| 273 |
plotInstructions.subPlotNo=2; |
|
| 274 |
plotInstructions.title=... |
|
| 275 |
['IC raster plot']; |
|
| 276 |
plotInstructions.figureNo=91; |
|
| 277 |
plotInstructions.displaydt=dtSpikes; |
|
| 278 |
plotInstructions.title='Brainstem 2nd level'; |
|
| 279 |
plotInstructions.yLabel='BF'; |
|
| 280 |
plotInstructions.yValues= savedBFlist; |
|
| 281 |
|
|
| 282 |
if size(ICoutput,1)>1 |
|
| 283 |
if sum(sum(ICoutput))<100 |
|
| 284 |
plotInstructions.rasterDotSize=3; |
|
| 285 |
end |
|
| 286 |
UTIL_plotMatrix(ICoutput, plotInstructions); |
|
| 287 |
end |
|
| 288 |
|
|
| 289 |
% plot signal (1) |
|
| 290 |
plotInstructions.displaydt=dt; |
|
| 291 |
plotInstructions.subPlotNo=1; |
|
| 292 |
plotInstructions.title=... |
|
| 293 |
['stimulus (Pa). ' num2str(signalRMSdb, '%4.0f') ' dB SPL']; |
|
| 294 |
UTIL_plotMatrix(savedInputSignal, plotInstructions); |
|
| 295 |
|
|
| 296 |
end |
|
| 250 | 297 |
|
| 251 | 298 |
%% figure(98) plot efferent control values as dB |
| 252 | 299 |
if showMapOptions.showEfferent |
| ... | ... | |
| 259 | 306 |
plotInstructions.zValuesRange= [-1 1]; |
| 260 | 307 |
plotInstructions.title= ['RMS level='... |
| 261 | 308 |
num2str(signalRMSdb, '%4.0f') ' dB SPL']; |
| 262 |
UTIL_plotMatrix(savedInputSignal', plotInstructions);
|
|
| 263 |
|
|
| 264 |
|
|
| 309 |
UTIL_plotMatrix(savedInputSignal, plotInstructions); |
|
| 310 |
|
|
| 311 |
|
|
| 265 | 312 |
plotInstructions.subPlotNo=2; |
| 266 | 313 |
plotInstructions.zValuesRange=[ -25 0]; |
| 267 | 314 |
plotInstructions.title= ['AR stapes attenuation (dB); tau='... |
| ... | ... | |
| 275 | 322 |
plotInstructions.yLabel= 'dB'; |
| 276 | 323 |
if strcmp(saveAN_spikesOrProbability,'spikes') |
| 277 | 324 |
rate2atten=DRNLParams.rateToAttenuationFactor; |
| 278 |
plotInstructions.title= ['MOC atten; tau=' ... |
|
| 279 |
num2str(DRNLParams.MOCtau) '; factor='... |
|
| 280 |
num2str(rate2atten, '%6.4f')]; |
|
| 325 |
plotInstructions.title= ['MOC atten; tau=' ...
|
|
| 326 |
num2str(DRNLParams.MOCtau) '; factor='...
|
|
| 327 |
num2str(rate2atten, '%6.4f')];
|
|
| 281 | 328 |
else |
| 282 | 329 |
rate2atten=DRNLParams.rateToAttenuationFactorProb; |
| 283 |
plotInstructions.title= ['MOC atten; tauProb=' ... |
|
| 284 |
num2str(DRNLParams.MOCtauProb) '; factor='... |
|
| 285 |
num2str(rate2atten, '%6.4f')]; |
|
| 330 |
plotInstructions.title= ['MOC atten; tauProb=' ...
|
|
| 331 |
num2str(DRNLParams.MOCtauProb) '; factor='...
|
|
| 332 |
num2str(rate2atten, '%6.4f')];
|
|
| 286 | 333 |
end |
| 287 | 334 |
plotInstructions.zValuesRange=[0 -DRNLParams.minMOCattenuationdB+5]; |
| 288 | 335 |
UTIL_plotMatrix(-20*log10(MOCattenuation), plotInstructions); |
| 289 | 336 |
hold on |
| 290 | 337 |
[r c]=size(MOCattenuation); |
| 291 |
if r>2 |
|
| 292 |
colorbar('southOutside')
|
|
| 338 |
if r>2 && showMapOptions.colorbar
|
|
| 339 |
colorbar('southOutside')
|
|
| 293 | 340 |
end |
| 294 | 341 |
set(plotInstructions.figureNo, 'name', 'AR/ MOC') |
| 295 |
|
|
| 342 |
|
|
| 296 | 343 |
binwidth=0.1; |
| 297 | 344 |
[PSTH ]=UTIL_PSTHmaker(20*log10(MOCattenuation), dt, binwidth); |
| 298 | 345 |
PSTH=PSTH*length(PSTH)/length(MOCattenuation); |
| 299 | 346 |
t=binwidth:binwidth:binwidth*length(PSTH); |
| 300 | 347 |
fprintf('\n\n')
|
| 301 |
% UTIL_printTabTable([t' PSTH']) |
|
| 302 |
% fprintf('\n\n')
|
|
| 303 |
|
|
| 348 |
% UTIL_printTabTable([t' PSTH'])
|
|
| 349 |
% fprintf('\n\n')
|
|
| 350 |
|
|
| 304 | 351 |
end |
| 305 | 352 |
|
| 306 |
%% ACF plot if required
|
|
| 353 |
%% ACF plot |
|
| 307 | 354 |
if showMapOptions.showACF |
| 308 | 355 |
tic |
| 309 |
method.dt=dt; |
|
| 310 |
method.segmentNo=1; |
|
| 311 |
method.nonlinCF=savedBFlist; |
|
| 312 |
|
|
| 313 |
minPitch= 80; maxPitch= 4000; numPitches=100; % specify lags |
|
| 314 |
pitches=10.^ linspace(log10(minPitch), log10(maxPitch),numPitches); |
|
| 315 |
pitches=fliplr(pitches); |
|
| 316 |
filteredSACFParams.lags=1./pitches; % autocorrelation lags vector |
|
| 317 |
filteredSACFParams.acfTau= .003; % time constant of running ACF |
|
| 318 |
filteredSACFParams.lambda= 0.12; % slower filter to smooth ACF |
|
| 319 |
filteredSACFParams.lambda= 0.01; % slower filter to smooth ACF |
|
| 320 |
|
|
| 321 |
filteredSACFParams.plotACFs=0; % special plot (see code) |
|
| 322 |
filteredSACFParams.plotFilteredSACF=0; % 0 plots unfiltered ACFs |
|
| 323 |
filteredSACFParams.plotMoviePauses=.3; % special plot (see code) |
|
| 324 |
|
|
| 325 |
filteredSACFParams.usePressnitzer=0; % attenuates ACF at long lags |
|
| 326 |
filteredSACFParams.lagsProcedure= 'useAllLags'; |
|
| 327 |
% filteredSACFParams.lagsProcedure= 'useBernsteinLagWeights'; |
|
| 328 |
% filteredSACFParams.lagsProcedure= 'omitShortLags'; |
|
| 329 |
filteredSACFParams.criterionForOmittingLags=3; |
|
| 330 |
filteredSACFParams.plotACFsInterval=200; |
|
| 331 |
|
|
| 332 | 356 |
if filteredSACFParams.plotACFs |
| 333 | 357 |
% plot original waveform on ACF plot |
| 334 | 358 |
figure(13), clf |
| ... | ... | |
| 339 | 363 |
title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']); |
| 340 | 364 |
end |
| 341 | 365 |
|
| 366 |
% compute ACF |
|
| 367 |
switch saveAN_spikesOrProbability |
|
| 368 |
case 'probability' |
|
| 369 |
inputToACF=ANprobRateOutput(end-length(savedBFlist)+1:end,:); |
|
| 370 |
otherwise |
|
| 371 |
inputToACF=ANoutput; |
|
| 372 |
end |
|
| 373 |
|
|
| 374 |
disp ('computing ACF...')
|
|
| 375 |
[P, BFlist, sacf, boundaryValue] = ... |
|
| 376 |
filteredSACF(inputToACF, dt, savedBFlist, filteredSACFParams); |
|
| 377 |
disp(' ACF done.')
|
|
| 378 |
|
|
| 342 | 379 |
% plot original waveform on summary/smoothed ACF plot |
| 343 | 380 |
figure(96), clf |
| 344 | 381 |
subplot(2,1,1) |
| ... | ... | |
| 347 | 384 |
xlim([0 t(end)]) |
| 348 | 385 |
title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']); |
| 349 | 386 |
|
| 350 |
|
|
| 351 |
% compute ACF |
|
| 352 |
switch saveAN_spikesOrProbability |
|
| 353 |
case 'probability' |
|
| 354 |
inputToACF=ANprobRateOutput.^0.5; |
|
| 355 |
otherwise |
|
| 356 |
inputToACF=ANoutput; |
|
| 357 |
end |
|
| 358 |
|
|
| 359 |
disp ('computing ACF...')
|
|
| 360 |
[P, BFlist, sacf, boundaryValue] = ... |
|
| 361 |
filteredSACF(inputToACF, method, filteredSACFParams); |
|
| 362 |
disp(' ACF done.')
|
|
| 363 |
|
|
| 364 |
% SACF |
|
| 387 |
% plot SACF |
|
| 388 |
figure(96) |
|
| 365 | 389 |
subplot(2,1,2) |
| 366 | 390 |
imagesc(P) |
| 391 |
% surf(filteredSACFParams.lags, t, P) |
|
| 367 | 392 |
ylabel('periodicities (Hz)')
|
| 368 | 393 |
xlabel('time (s)')
|
| 369 |
title(['running smoothed (root) SACF. ' saveAN_spikesOrProbability ' input'])
|
|
| 370 |
pt=[1 get(gca,'ytick')]; % force top xtick to show
|
|
| 394 |
title(['running smoothed SACF. ' saveAN_spikesOrProbability ' input']) |
|
| 395 |
pt=[1 get(gca,'ytick')]; % force top ytick to show
|
|
| 371 | 396 |
set(gca,'ytick',pt) |
| 397 |
pitches=1./filteredSACFParams.lags; % autocorrelation lags vector |
|
| 372 | 398 |
set(gca,'ytickLabel', round(pitches(pt))) |
| 399 |
[nCH nTimes]=size(P); |
|
| 400 |
t=dt:dt:dt*nTimes; |
|
| 373 | 401 |
tt=get(gca,'xtick'); |
| 374 | 402 |
set(gca,'xtickLabel', round(100*t(tt))/100) |
| 375 | 403 |
end |
| ... | ... | |
| 395 | 423 |
% end |
| 396 | 424 |
|
| 397 | 425 |
function ANsmooth = makeANsmooth(ANresponse, sampleRate, winSize, hopSize) |
| 398 |
if nargin < 3
|
|
| 399 |
winSize = 25; %default 25 ms window
|
|
| 400 |
end
|
|
| 401 |
if nargin < 4
|
|
| 402 |
hopSize = 10; %default 10 ms jump between windows
|
|
| 403 |
end
|
|
| 404 |
|
|
| 405 |
winSizeSamples = round(winSize*sampleRate/1000);
|
|
| 406 |
hopSizeSamples = round(hopSize*sampleRate/1000);
|
|
| 407 |
|
|
| 408 |
% smooth
|
|
| 409 |
hann = hanning(winSizeSamples);
|
|
| 410 |
|
|
| 411 |
ANsmooth = [];%Cannot pre-allocate a size as it is unknown until the enframing
|
|
| 412 |
for chan = 1:size(ANresponse,1)
|
|
| 413 |
f = enframe(ANresponse(chan,:), hann, hopSizeSamples);
|
|
| 414 |
ANsmooth(chan,:) = mean(f,2)'; %#ok<AGROW> see above comment
|
|
| 415 |
end
|
|
| 426 |
if nargin < 3 |
|
| 427 |
winSize = 25; %default 25 ms window |
|
| 428 |
end |
|
| 429 |
if nargin < 4 |
|
| 430 |
hopSize = 10; %default 10 ms jump between windows |
|
| 431 |
end |
|
| 432 |
|
|
| 433 |
winSizeSamples = round(winSize*sampleRate/1000); |
|
| 434 |
hopSizeSamples = round(hopSize*sampleRate/1000); |
|
| 435 |
|
|
| 436 |
% smooth |
|
| 437 |
hann = hanning(winSizeSamples); |
|
| 438 |
|
|
| 439 |
ANsmooth = [];%Cannot pre-allocate a size as it is unknown until the enframing |
|
| 440 |
for chan = 1:size(ANresponse,1) |
|
| 441 |
f = enframe(ANresponse(chan,:), hann, hopSizeSamples); |
|
| 442 |
ANsmooth(chan,:) = mean(f,2)'; %#ok<AGROW> see above comment |
|
| 443 |
end |
|
| 416 | 444 |
% end% ------ OF makeANsmooth |
| utilities/stimulusCreate.m | ||
|---|---|---|
| 399 | 399 |
stimulus=applyGaussianRamps(stimulus, -stim.rampOnDur, 1/dt); |
| 400 | 400 |
end |
| 401 | 401 |
|
| 402 |
% Initial silence
|
|
| 402 |
% begin silence
|
|
| 403 | 403 |
% start with a signal of the right length consisting of zeros |
| 404 | 404 |
signal=zeros(1, globalStimParams.nSignalPoints); |
| 405 | 405 |
% identify start of stimulus |
| ... | ... | |
| 414 | 414 |
|
| 415 | 415 |
% time signature |
| 416 | 416 |
time=dt:dt:dt*length(signal); |
| 417 |
% figure(22), plot(signal), title([num2str(ear) ' - ' num2str(componentNo)]),pause (1)
|
|
| 417 |
% figure(22), plot(signal), title([num2str(ear) ' - ' num2str(componentNo)]),pause (1) |
|
| 418 | 418 |
|
| 419 | 419 |
try |
| 420 | 420 |
% create a column vector and trap if no signal has been created |
| ... | ... | |
| 574 | 574 |
function stimulus=makeOHIOtones(globalStimParams, stimComponents) |
| 575 | 575 |
|
| 576 | 576 |
% Generates a stimulus consisting of one or more 10-ms tones |
| 577 |
% Tones are presented at 10-ms intervals |
|
| 577 |
% The length of the list of frequencies determines the numberof tones |
|
| 578 |
% Tones are either presented at 10-ms intervals or simultaneously |
|
| 579 |
% all tones are individually ramped |
|
| 578 | 580 |
% Each tone has its own amplitude and its own ramp. |
| 579 | 581 |
|
| 580 | 582 |
frequencies=stimComponents.frequencies; |
| 581 | 583 |
amplitudesdB=stimComponents.amplitudesdB; |
| 582 | 584 |
nFrequencies=length(frequencies); |
| 583 | 585 |
|
| 586 |
if amplitudesdB==-100 |
|
| 587 |
% catch trial |
|
| 588 |
amplitudesdB=repmat(-100,1,nFrequencies); |
|
| 589 |
end |
|
| 590 |
|
|
| 584 | 591 |
dt=globalStimParams.dt; |
| 592 |
|
|
| 585 | 593 |
toneDuration=.010; |
| 586 | 594 |
time=dt:dt:toneDuration; |
| 587 | 595 |
|
| 588 |
rampDuration=stimComponents.rampOnDur; |
|
| 596 |
% fixed ramp, silenceDuration, toneDuration |
|
| 597 |
rampDuration=0.005; |
|
| 589 | 598 |
rampTime=dt:dt:rampDuration; |
| 590 |
ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ones(1,length(time)-length(rampTime))]; |
|
| 599 |
ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ... |
|
| 600 |
ones(1,length(time)-length(rampTime))]; |
|
| 591 | 601 |
ramp=ramp.*fliplr(ramp); |
| 592 | 602 |
|
| 593 |
stimulus= zeros(1,round((toneDuration+globalStimParams.beginSilences(end))/dt)+1); |
|
| 603 |
silenceDuration=0.010; |
|
| 604 |
silenceDurationLength=round(silenceDuration/dt); |
|
| 605 |
initialSilence=zeros(1,silenceDurationLength); |
|
| 606 |
|
|
| 607 |
silenceToneDuration=toneDuration + silenceDuration; |
|
| 608 |
silenceToneDurationLength=round(silenceToneDuration/dt); |
|
| 609 |
|
|
| 610 |
% OHIO spect produces simultaneous tones |
|
| 611 |
if strcmp(stimComponents.OHIOtype,'OHIOspect') |
|
| 612 |
totalDuration=silenceToneDuration; |
|
| 613 |
else |
|
| 614 |
totalDuration=silenceToneDuration*nFrequencies; |
|
| 615 |
end |
|
| 616 |
|
|
| 617 |
totalDurationLength=round(totalDuration/dt); |
|
| 618 |
stimulus=zeros(1,totalDurationLength); |
|
| 619 |
toneBeginPTR=1; |
|
| 594 | 620 |
|
| 595 | 621 |
for i=1:nFrequencies |
| 596 |
toneBeginPTR=round(globalStimParams.beginSilences(i)/dt)+1; |
|
| 597 |
|
|
| 598 | 622 |
frequency=frequencies(i); |
| 599 | 623 |
dBSPL=amplitudesdB(i); |
| 600 | 624 |
amplitude=28e-6* 10.^(dBSPL/20); |
| 601 | 625 |
tone=amplitude*sin(2*pi*frequency*time); |
| 602 | 626 |
tone=tone.*ramp; |
| 603 |
stimulus(toneBeginPTR:toneBeginPTR+length(tone)-1)=stimulus(toneBeginPTR:toneBeginPTR+length(tone)-1)+tone; |
|
| 627 |
% stimulus is normally zeros except for OHIOspect |
|
| 628 |
stimulus(toneBeginPTR:toneBeginPTR+silenceToneDurationLength-1)=... |
|
| 629 |
[initialSilence tone]+... |
|
| 630 |
stimulus(toneBeginPTR:toneBeginPTR+silenceToneDurationLength-1); |
|
| 631 |
if ~strcmp(stimComponents.OHIOtype,'OHIOspect') |
|
| 632 |
toneBeginPTR=toneBeginPTR+silenceToneDurationLength; |
|
| 633 |
end |
|
| 604 | 634 |
end |
| 605 | 635 |
% figure(2), plot( stimulus') |
| 606 | 636 |
|
| utilities/temp.m | ||
|---|---|---|
| 1 |
frequencies=[1000 1250]; |
|
| 2 |
amplitudesdB=[20 23]; |
|
| 3 |
nFrequencies=length(frequencies); |
|
| 4 |
|
|
| 5 |
dt=0.0001; |
|
| 6 |
|
|
| 7 |
toneDuration=.010; |
|
| 8 |
time=dt:dt:toneDuration; |
|
| 9 |
|
|
| 10 |
% fixed ramp, silenceDuration, toneDuration |
|
| 11 |
rampDuration=0.005; |
|
| 12 |
rampTime=dt:dt:rampDuration; |
|
| 13 |
ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ... |
|
| 14 |
ones(1,length(time)-length(rampTime))]; |
|
| 15 |
ramp=ramp.*fliplr(ramp); |
|
| 16 |
|
|
| 17 |
silenceDuration=0.010; |
|
| 18 |
silenceDurationLength=round(silenceDuration/dt); |
|
| 19 |
initialSilence=zeros(1,silenceDurationLength); |
|
| 20 |
|
|
| 21 |
silenceToneDuration=toneDuration + silenceDuration; |
|
| 22 |
silenceToneDurationLength=round(silenceToneDuration/dt); |
|
| 23 |
|
|
| 24 |
totalDuration=silenceToneDuration*nFrequencies; |
|
| 25 |
totalDurationLength=round(totalDuration/dt); |
|
| 26 |
stimulus=zeros(1,totalDurationLength); |
|
| 27 |
toneBeginPTR=1; |
|
| 28 |
|
|
| 29 |
for i=1:nFrequencies |
|
| 30 |
frequency=frequencies(i); |
|
| 31 |
dBSPL=amplitudesdB(i); |
|
| 32 |
amplitude=28e-6* 10.^(dBSPL/20); |
|
| 33 |
tone=amplitude*sin(2*pi*frequency*time); |
|
| 34 |
tone=tone.*ramp; |
|
| 35 |
stimulus(toneBeginPTR:toneBeginPTR+silenceToneDurationLength-1)=... |
|
| 36 |
[initialSilence tone]; |
|
| 37 |
toneBeginPTR=toneBeginPTR+silenceToneDurationLength; |
|
| 38 |
end |
|
| 39 |
figure(2), plot( stimulus') |
|
Also available in: Unified diff