pulsetrain.m

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

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