rmeddis@38: function wholeNerveCAP  = UTIL_CAPgenerator...
rmeddis@38:     (ANresponse, dt, BFlist, numANfibers, plotCAP)
rmeddis@38: % 
rmeddis@38: % Generates a compound action potential by convolving an impulse repsonse,
rmeddis@38: % as defined by Mark Chertoff (2004, JASA), with the response of the 
rmeddis@38: % auditory nerve.
rmeddis@38: %
rmeddis@38: %
rmeddis@38: % -e(-k*time)*SIN(2*PI()*f*time)
rmeddis@38: %
rmeddis@38: % mu(t) = e^-kt * sin(omega*t)
rmeddis@38: % omega = 2 * pi * f
rmeddis@38: %
rmeddis@38: %
rmeddis@38: % Robert T. Ferry
rmeddis@38: % 01st May 2008
rmeddis@38: %
rmeddis@38: %
rmeddis@38: 
rmeddis@38: nChannels=length(BFlist);
rmeddis@38: [r nSpikeEpochs]=size(ANresponse);
rmeddis@38: 
rmeddis@38: wholeNerveCAP       = [];
rmeddis@38: channelCAP          = [];
rmeddis@38: e                   = exp(1);
rmeddis@38: k                   = 1000;
rmeddis@38: impulseDuration     = 0.01;
rmeddis@38: impulseFrequency    = 750;
rmeddis@38: impulseTime         = dt:dt:impulseDuration;
rmeddis@38: impulseResponse     = -e.^(-k*impulseTime).*sin(2*pi*impulseFrequency*impulseTime);
rmeddis@38: impulseResponse=impulseResponse-mean(impulseResponse);
rmeddis@38: 
rmeddis@38: % WholeNerveCAP
rmeddis@38: ANoutput = sum(ANresponse, 1);
rmeddis@38: convolvedWholeNerveCAP = conv(ANoutput, impulseResponse(1,:));
rmeddis@38: % truncate
rmeddis@38: convolvedWholeNerveCAP=convolvedWholeNerveCAP(1:nSpikeEpochs);
rmeddis@38: 
rmeddis@38: % apply measurement time constant
rmeddis@38: sampleRate=1/dt;
rmeddis@38: upperFreq=sampleRate/4;
rmeddis@38: lowPassCutOff=40;
rmeddis@38: wholeNerveCAP = UTIL_Butterworth(convolvedWholeNerveCAP, dt, lowPassCutOff, upperFreq, 1);
rmeddis@38: % or do not do this
rmeddis@38: % wholeNerveCAP = convolvedWholeNerveCAP;
rmeddis@38: 
rmeddis@38: % Plot output?
rmeddis@38: 
rmeddis@38: CAPtime = dt:dt:dt*length(wholeNerveCAP);
rmeddis@38: 
rmeddis@38: if (plotCAP == 1)
rmeddis@38:     figure(9)
rmeddis@38: 
rmeddis@38:     subplot(3,1,1), plot(impulseTime, impulseResponse)
rmeddis@38:     title('Impulse response')
rmeddis@38:     xlabel('Time (s)'), ylabel('Amplitude (Pa)')
rmeddis@38:     xlim([0 max(impulseTime)]), ylim([-inf inf])
rmeddis@38: 	
rmeddis@38:     subplot(3,1,2), plot(CAPtime, wholeNerveCAP)
rmeddis@38:     title(['AN CAP (whole-nerve)  '   num2str(length(BFlist)) ' channels' num2str(numANfibers) ' fibers/ch'])
rmeddis@38:     xlabel('Time (s)'), ylabel('Amplitude (Pa)')
rmeddis@38:     xlim([0 max(CAPtime)]), ylim([-inf inf])   
rmeddis@38: 	
rmeddis@38: end