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