www.gusucode.com > 时间序列工具箱 > 时间序列工具箱/timeseries_xg/matlab_download/xg_timeseries/cochlear/pulsetrain.m
function [prperpulse, stdperpulse]=pulsetrain(stim,properties,pulserate,pulsewidth) % [prperpulse, stdperpulse]=pulsetrain(stim,properties,pulserate,pulsewidth) % % INPUT: % % stim is the stimulus intensity in units of microA, NOT dB % % properties is a structure containing the model parameters: % % properties.resptype is the fiber I/O function type: % 'stp' = deterministic (step-function) model % '3_p' = 3-piece linear model % 'erf' = stochastic (error-function) model % % properties.thrsh is the matrix of fiber thresholds % % properties.noisestd is the matrix of fiber noise standard deviations % % pulserate is the rate of stimulation in pulses/s % % pulsewidth is the phase duration in microseconds/phase % % OUTPUT: % % prperpulse is the discharge probability (mean discharge rate) per pulse % % stdperpulse is the standard deviation in discharge rate per pulse % % Usage Agreement: Any publication of results obtained using this model should cite both- % % [1] Bruce, I. C., Irlicht, L. S., White, M. W., O'Leary, S. J., Dynes, % S., Javel, E., Clark, G. M. (1999) "A stochastic model of the % electrically stimulated auditory nerve: Pulse-train response," % IEEE Transactions on Biomedical Engineering, 46(6):630-637. % % and % % [2] Bruce, I. C. (1997) Spatiotemporal coding of sound in the auditory % nerve for cochlear implants, PhD Thesis, The University of Melbourne, Australia. % Copyright (c) 1996-1999 Ian Bruce % Created by Ian Bruce, 1996 % Released for public distribution as version 1.1, 1999 Fs = 1e5; % Sampling frequency in Hz resptype = properties.resptype; thrsh = properties.thrsh; noisestd = properties.noisestd; % Run pulse-train model for each individual fiber at each stimulus intensity for i=1:size(stim,1) for j=1:size(stim,2) disp([num2str(pulserate) 'Hz - Run ' num2str(size(stim,2)*(i-1)+j) '/' num2str(size(stim,1)*size(stim,2))]) [prperpulse(i,j), stdperpulse(i,j)] = pulsetrain_perelement(stim(i,j),resptype,thrsh(i,j),noisestd(i,j),pulserate,pulsewidth,Fs); end end function [prperpulse, stdperpulse, isiprs, Pinf]=pulsetrain_perelement(level,resptype,thrsh,noisestd,pulserate,pulsewidth,Fs) % % [prperpulse, stdperpulse, isiprs, Pinf]=pulsetrain_perelement(level,resptype,thrsh,noisestd,pulserate,pulsewidth,Fs) % For Vrefr relative to Vthr absref = 0.7e-3; % Absolute refractory period in seconds if pulserate>1/absref error(['this model is not appropriate for pulserates above ' num2str(round(1/absref)) ' pulses/s']) end relref = 20e-3; % Relative refractory period in seconds magn = 0.97; % Relative magnitude of refactory function tconst = 1.32e-3; % Time-constant of refactory function in seconds y0 = 1; % Asymptotic relative magnitude of refactory function dt = 1/Fs; % bin size in seconds pulsewidth = pulsewidth*1e-6; % convert from microseconds/phase to seconds/phase % Need to calculate number of bins per pulse, to calculate mean and variance of renewal time % depending on which bin the previous spike was in, i.e., the 'starting' bin starts=1:round(pulsewidth*Fs); % Number of starting bins stlst = length(starts); % Initialize 3D refractory function matrix r = ones(round(relref/dt),1,stlst); % Calculate refractory function for each starting bin for lp=1:stlst r(1:round(absref/dt)-1+starts(lp),1,starts(lp))=Inf; r(round(absref/dt)+starts(lp):round(relref/dt),1,starts(lp))=y0+magn*exp(-([round(absref/dt)+starts(lp):round(relref/dt)]-(round(absref/dt)+starts(lp)-1))/tconst*dt); % exp r end % Calculate number of bins per stimulus period bins = round(1/pulserate/dt); % Pad out refractory function to integer multiple of the stimulus period rtmp = cat(1,r,ones(bins*round(length(r)/bins+1)-length(r),1,stlst)); rlen = size(rtmp,1); % Reshape refractory function matrix so that: % 1st dimension is the bin number within a stimulus cycle (from the start of one pulse to the start of the next), % 2nd dimension is the stimulus cycle, and % 3rd dimension is the starting bin. rftemp = reshape(rtmp,bins,rlen/bins,stlst); % Get refractory function for bins that fall within a pulse from the 2nd pulse onwards % (the 'previous' spike has occurred in the 1st pulse) rf = rftemp(1:round(pulsewidth*Fs),2:rlen/bins,starts); dims = size(rf); % Number of pulses, excluding the 1st pulse where the 'previous' spike has occurred plss = dims(2); % Initialize matrices f = zeros(plss,1,stlst); p = zeros(dims); stim = ones(dims); % Put model parameters into a single structure properties.resptype = resptype; properties.thrsh = thrsh; properties.noisestd = noisestd; % Calculate asymptotic discharge probability using single-pulse model Pinf = singlepulse(level,properties); % Create matrix of stimulus level levels=stim.*level; % Create matrix of refractory-modified threshold properties.thrsh = thrsh.*rf; % Calculate discharge probabilities using single-pulse model p = singlepulse(levels,properties); % Get discharge probability for each pulse P = p(dims(1),:,:); % Transpose the 3D matrix P = permute(P,[2 1 3]); q = 1-P; Q = cumprod(q); % See Lemma 7.3.2 of [2] f(1,1,:) = P(1,1,:); % Eqn. (7.7) of [2] f(2:plss,1,:) = P(2:plss,1,:).*Q(1:plss-1,1,:); % Eqn. (7.7) of [2] % Calculate pulse numbers ks=repmat([1:plss]',[1 1 stlst]); if Pinf==0 % Handle case where asymptotic discharge probability is zero tmean = Inf; disrate = 0; tvar = 0; ratestd = 0; isipr = zeros(size(f)); else tmean = sum(ks.*f)+Q(end,1,:).*(plss+1./Pinf); % Eqn. (7.9) of [2] disrate = 1./tmean; % Eqn. (7.17) of [2] tvar = sum(ks.^2.*f)+Q(end,1,:).*(-2*(plss+1)-3./Pinf+1+(plss+1)^2+2*(plss+1)./Pinf+2./Pinf.^2)-tmean.^2; % Eqn. (7.10) of [2] ratestd = sqrt(tvar./(tmean).^3); % Eqn. (7.18) of [2] isipr = f; end % Calculate matrix B using Eqn. (7.11) of [2] X=diff(cat(1,zeros(1,plss,stlst),p)); % Eqn. (7.12) of [2] Y=permute(repmat(cat(1,ones(1,1,stlst),Q(1:plss-1,:,:)),[1 dims(1) 1]),[2 1 3]); Z=Q(end,1,:); AA=cat(1, ones(1,1,stlst), zeros(dims(1)-1,1,stlst)); Barray = sum((X.*Y),2)+arraydotprod(AA,Z); % Reduce to 2D matrix B = squeeze(Barray); [v, d] = eig(B); % Find eigenvectors and eigenvalues [i, j] = find(round(d)==1); % Find which eigenvector has eigenvalue 1 b = v(:,j)/sum(v(:,j)); % Normalize eigenvector % Turn vectors into column vectors disrate = disrate(:); ratestd = ratestd(:); prperpulse = sum(disrate.*b); % Eqn. (7.15) of [2] stdperpulse = sqrt(sum(ratestd.^2.*b)); % Eqn. (7.16) of [2] isiprs = squeeze(isipr)*b; % Eqn. (7.19) of [2] function Y=arraydotprod(A,B); if ndims(A)~=3 Y=[]; disp('??? Error using ==> arraydotprod') disp('First input must be 3-D array') elseif ndims(A)~=3 Y=[]; disp('??? Error using ==> arraydotprod') disp('Second input must be 3-D array') elseif size(A,2)~=size(B,1) Y=[]; disp('??? Error using ==> arraydotprod') disp('Inner matrix dimensions must agree.') elseif size(A,3)~=size(B,3) Y=[]; disp('??? Error using ==> arraydotprod') disp('3rd dimensions must agree.') elseif size(B,2)~=1 Y=[]; disp('??? Error using ==> arraydotprod') disp('Second input must have singleton 2nd dimension.') else Y=sum(A.*repmat(permute(B,[2 1 3]),[size(A,1) 1 1]),2); end