Daniel@0: function y = MeddisHairCell(data,sampleRate,subtractSpont)
Daniel@0: % y = MeddisHairCell(data,sampleRate)
Daniel@0: % This function calculates Ray Meddis' hair cell model for a
Daniel@0: % number of channels.  Data is arrayed as one channel per row.
Daniel@0: % All channels are done in parallel (but each time step is
Daniel@0: % sequential) so it will be much more efficient to process lots
Daniel@0: % of channels at once.
Daniel@0: 
Daniel@0: % (c) 1998 Interval Research Corporation  
Daniel@0: % Changed h and added comment at suggestion of Alain de Cheveigne'   12/11/98
Daniel@0: 
Daniel@0: if (nargin<3),  subtractSpont=0;  end
Daniel@0: 
Daniel@0: % Parameters from Meddis' April 1990 JASA paper.
Daniel@0: M = 1;
Daniel@0: A = 5;
Daniel@0: B = 300;
Daniel@0: g = 2000;
Daniel@0: y = 5.05;
Daniel@0: l = 2500;
Daniel@0: r = 6580;
Daniel@0: x = 66.31;
Daniel@0: h = 50000;	% This parameter scales the discharge rate. Adjust as necessary.
Daniel@0: 		% In combination with the gammatone filterbank (ERBFilterBank),
Daniel@0: 		% h=50000 will produce a steady-state average discharge
Daniel@0: 		% probability of about 135 spikes/s within the 1kHz channel,
Daniel@0: 		% for an input consisting of a 1 kHz sinewave at 60 dB SPL
Daniel@0: 		% (0 dB SPL corresponds to an RMS level of 1.0 at the
Daniel@0: 		% input of the gammatone filter).  Scaling and constant 
Daniel@0: 		% courtesy of Alain de Cheveigne'
Daniel@0: 
Daniel@0: 
Daniel@0: % Internal constants
Daniel@0: dt = 1/sampleRate;
Daniel@0: gdt = g*dt;
Daniel@0: ydt = y*dt;
Daniel@0: ldt = l*dt;
Daniel@0: rdt = r*dt;
Daniel@0: xdt = x*dt;
Daniel@0: [numChannels dataLength] = size(data);
Daniel@0: 
Daniel@0: % Initial values 
Daniel@0: kt = g*A/(A+B);
Daniel@0: spont = M*y*kt/(l*kt+y*(l+r));
Daniel@0: c = spont * ones(numChannels,1);
Daniel@0: q = c*(l+r)/kt;
Daniel@0: w = c*r/x;
Daniel@0: zeroVector = zeros(numChannels,1);
Daniel@0: 
Daniel@0: % Now iterate through each time slice of the data.  Use the
Daniel@0: % max function to implement the "if (0>" test.
Daniel@0: y = zeros(numChannels, dataLength);
Daniel@0: for i = 1:dataLength
Daniel@0:  limitedSt = max(data(:,i)+A,0);
Daniel@0:  kt = gdt*limitedSt./(limitedSt+B);
Daniel@0:  replenish = max(ydt*(M-q),zeroVector);
Daniel@0:  eject = kt.*q;
Daniel@0:  loss = ldt.*c;
Daniel@0:  reuptake = rdt.*c;
Daniel@0:  reprocess = xdt.*w;
Daniel@0: 
Daniel@0:  q = q + replenish - eject + reprocess;
Daniel@0:  c = c + eject - loss - reuptake;
Daniel@0:  w = w + reuptake - reprocess;
Daniel@0:  y(:,i) = c;
Daniel@0: end
Daniel@0: 
Daniel@0: y = h .* y;
Daniel@0: 
Daniel@0: if (subtractSpont > 0)
Daniel@0:  y=max(0,y-spont);
Daniel@0: end