totalan.m

这是时间序列的matlab程序应用

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