www.gusucode.com > 时间序列工具箱 > 时间序列工具箱/timeseries_xg/matlab_download/xg_timeseries/cochlear/totalan.m
function [cnt, cntstd, totalancnt, totalanstd]=totalan(stim,numofpulses,pulserate,pulsewidth,resptype,electconfig,N) % [cnt, cntstd, totalancnt, totalanstd]=totalan(stim,numofpulses,pulserate,pulsewidth,resptype,electconfig,N) % % INPUT: % % stim is a vector of stimulus intensities in dB re. 1 microA % % numofpulses is the number of pulses in the pulse train % % pulserate is the rate of stimulation in pulses/s % note that for pulserate <= 100 pulse/s the singlepulse model is used, % otherwise the pulsetrain model is used % % pulsewidth is the phase duration in microseconds/phase % % resptype is the fiber I/O function type: % 'stp' = deterministic (step-function) model % '3_p' = 3-piece linear model % 'erf' = stochastic (error-function) model % % electconfig is the electrode configuration: % 'mp' = monopolar; 'bp' = bipolar % % N is the number of fibers % % OUTPUT: % % cnt is the mean number of discharges for each fiber % cntstd is the standard deviation in number of discharges for each fiber % % totalcnt is the mean number of discharges for the total Auditory Nerve (AN) % totalcntstd is the standard deviation in number of discharges for the total AN % % e.g., [cnt, cntstd, totalancnt, totalanstd]=totalan(40:2:70,1,1,100,'erf','mp',1e4) % gives the model output for stimulus intensities from 40 dB to 70 dB re. 1 microA % in 2 dB steps, for a single pulse of phase duration 100 microseconds/phase, % the stochastic model, monopolar stimulation and 10,000 fibers. % % e.g., [cnt, cntstd, totalancnt, totalanstd]=totalan(45:5:65,10,400,200,'stp','bp',1e3) % gives the model output for stimulus intensities from 45 dB to 55 dB re. 1 microA % in 1 dB steps, for 10 pulses at 400 pulses/s and phase duration 200 microseconds/phase, % the deterministic model, bipolar stimulation and 1,000 fibers. % % Usage Agreement: Any publication of results obtained using this model should cite- % % [1] Bruce, I. C., White, M. W., Irlicht, L. S., O'Leary, S. J., Dynes, % S., Javel, E., Clark, G. M. (1999) "A stochastic model of the % electrically stimulated auditory nerve: Single-pulse response," % IEEE Transactions on Biomedical Engineering, 46(6):617-629. % Copyright (c) 1996-1999 Ian Bruce % Created by Ian Bruce, 1996 % Released for public distribution as version 1.1, 1999 % Check for valid stimulus vector if ndims(stim)>2 error('stim has >2 dimensions - it should be a vector') elseif (size(stim,1)~=1)&(size(stim,2)~=1) error('stim is a matrix - it should be a vector') elseif (size(stim,1)~=1)&(size(stim,2)==1) stim = stim'; end % Check for valid number of pulses if ndims(numofpulses)>2 error('numofpulses has >2 dimensions - it should be a scalar') elseif (size(numofpulses,1)~=1)|(size(numofpulses,2)~=1) error('numofpulses is a matrix (or vector) - it should be a scalar') elseif numofpulses<1 error('numofpulses is <1 - it should be non-zero and positive') end % Check for valid pulse rate if ndims(pulserate)>2 error('pulserate has >2 dimensions - it should be a scalar') elseif (size(pulserate,1)~=1)|(size(pulserate,2)~=1) error('pulserate is a matrix (or vector) - it should be a scalar') elseif pulserate<1 error('pulserate is <1 - it should be non-zero and positive') end % Check for valid phase duration if ndims(pulsewidth)>2 error('pulsewidth has >2 dimensions - it should be a scalar') elseif (size(pulsewidth,1)~=1)|(size(pulsewidth,2)~=1) error('pulsewidth is a matrix (or vector) - it should be a scalar') elseif pulsewidth<1 error('pulsewidth is <1 - it should be non-zero and positive') end % Check for valid electrode configuration if electconfig == 'mp' % Monopolar electrode configuration rateatten = 15; % Attenuation of 0.5 dB/mm for a 30 mm cochlea elseif electconfig == 'bp' % Bipolar electrode configuration rateatten = 120; % Attenuation of 4.0 dB/mm for a 30 mm cochlea else error(['electconfig ' electconfig ' not supported']) end % Check for valid number of fibers if ndims(N)>2 error('N has >2 dimensions - it should be a scalar') elseif (size(N,1)~=1)|(size(N,2)~=1) error('N is a matrix (or vector) - it should be a scalar') elseif N<1 error('N is <1 - it should be non-zero and positive') end % Put stim through current spread model if N == 1 stimmat = stim; % No attenuation for a single fiber, i.e., assume it is the fiber closest to the electrode else stimmat = ispread(stim,0.5,N,rateatten); % electrode is 0.5 of the way into the cochlea % ispread function given below end % Calculate threshold for each fiber thrsh=121.0354*pulsewidth^(-0.1821); % Eqn. (5) of [1]; threshold is at Pr of 0.5 thresholdrange = 10; % dB rand('state',0); % Initialize seed of random variable for thresholdshifts thresholdshifts = thresholdrange.*rand(N,1)*ones(1,length(stim))-thresholdrange/2; % dB % Calculate Relative Spread (RS) for each fiber rsmean = 0.1222 + 9.5063e-5*pulsewidth - 7.9042e-9*pulsewidth^2; % Eqn. (6) of [1] rsstd = 0.06; randn('state',pi); % Initialize seed of random variable for relativespread relativespread=rsstd.*randn(N,1)*ones(1,length(stim))+rsmean; % RS of 0 can give divide by zero error in stochastic model, so make 0s into smallest possible floating point relativespread=(relativespread>0).*relativespread+(relativespread<=0).*eps; % Convert RS to noise standard deviations; note that standard deviations needs to be given in units of microA noisestd=relativespread.*10.^((thrsh+thresholdshifts)/20); % Eqn. (4) of [1] % Initialize response matrices response = zeros(size(stimmat)); responsestd = zeros(size(stimmat)); % Put model parameters into a single structure; note that thresholds need to be given in units of microA properties = struct('resptype',resptype,'thrsh',10.^((thrsh+thresholdshifts)/20),'noisestd',noisestd); % Run single-fiber model; note that stimulus needs to be in units of microA if pulserate<=100 % Use single-pulse model response = singlepulse(10.^(stimmat/20),properties); responsestd = sqrt(response.*(1-response)); else % Use pulse-train model [response, responsestd] = pulsetrain(10.^(stimmat/20),properties,pulserate,pulsewidth); end % Single-fiber models give output as discharge mean and standard deviation per pulse, % so we need to take into account number of pulses here cnt = numofpulses*response; cntstd = sqrt(numofpulses)*responsestd; % Calculate summed responses totalancnt = sum(cnt,1); totalanstd = sqrt(sum(cntstd.^2,1)); function Ivect=ispread(IdB,site,N,rateatten) % IdB is the stimulus current vector % site is the electrode placement within the cochlea % N is the number of fibers % rateatten is the attenuation rate in dB across the whole cochlea % Check for valid electrode placement if ndims(site)>2 error('site has >2 dimensions - it should be a scalar') elseif (size(site,1)~=1)|(size(site,2)~=1) error('site is a matrix (or vector) - it should be a scalar') elseif (site<0)|(site>1) error('site is <0 or >1 - it should >=0 and <=1') end % Check for valid attenuation rate if ndims(rateatten)>2 error('rateatten has >2 dimensions - it should be a scalar') elseif (size(rateatten,1)~=1)|(size(rateatten,2)~=1) error('rateatten is a matrix (or vector) - it should be a scalar') elseif rateatten<=0 error('rateatten is <=0 - it should be positive') end Ivect = ones(N,1)*IdB; % Replicate Idb for each fiber atten = zeros(N,1)*ones(size(IdB)); % Initialize attenuation matrix % Calulculate attenuation for each fiber atten = atten + abs(rateatten*(1/(N-1)*[0:N-1]'-site))*ones(size(IdB)); Ivect = Ivect - atten; % Attenuate the stimulus vector