comparison toolboxes/MIRtoolbox1.3.2/AuditoryToolbox/MeddisHairCell.m @ 0:e9a9cd732c1e tip

first hg version after svn

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