|
wolffd@0
|
1 function y = MeddisHairCell(data,sampleRate,subtractSpont)
|
|
wolffd@0
|
2 % y = MeddisHairCell(data,sampleRate)
|
|
wolffd@0
|
3 % This function calculates Ray Meddis' hair cell model for a
|
|
wolffd@0
|
4 % number of channels. Data is arrayed as one channel per row.
|
|
wolffd@0
|
5 % All channels are done in parallel (but each time step is
|
|
wolffd@0
|
6 % sequential) so it will be much more efficient to process lots
|
|
wolffd@0
|
7 % of channels at once.
|
|
wolffd@0
|
8
|
|
wolffd@0
|
9 % (c) 1998 Interval Research Corporation
|
|
wolffd@0
|
10 % Changed h and added comment at suggestion of Alain de Cheveigne' 12/11/98
|
|
wolffd@0
|
11
|
|
wolffd@0
|
12 if (nargin<3), subtractSpont=0; end
|
|
wolffd@0
|
13
|
|
wolffd@0
|
14 % Parameters from Meddis' April 1990 JASA paper.
|
|
wolffd@0
|
15 M = 1;
|
|
wolffd@0
|
16 A = 5;
|
|
wolffd@0
|
17 B = 300;
|
|
wolffd@0
|
18 g = 2000;
|
|
wolffd@0
|
19 y = 5.05;
|
|
wolffd@0
|
20 l = 2500;
|
|
wolffd@0
|
21 r = 6580;
|
|
wolffd@0
|
22 x = 66.31;
|
|
wolffd@0
|
23 h = 50000; % This parameter scales the discharge rate. Adjust as necessary.
|
|
wolffd@0
|
24 % In combination with the gammatone filterbank (ERBFilterBank),
|
|
wolffd@0
|
25 % h=50000 will produce a steady-state average discharge
|
|
wolffd@0
|
26 % probability of about 135 spikes/s within the 1kHz channel,
|
|
wolffd@0
|
27 % for an input consisting of a 1 kHz sinewave at 60 dB SPL
|
|
wolffd@0
|
28 % (0 dB SPL corresponds to an RMS level of 1.0 at the
|
|
wolffd@0
|
29 % input of the gammatone filter). Scaling and constant
|
|
wolffd@0
|
30 % courtesy of Alain de Cheveigne'
|
|
wolffd@0
|
31
|
|
wolffd@0
|
32
|
|
wolffd@0
|
33 % Internal constants
|
|
wolffd@0
|
34 dt = 1/sampleRate;
|
|
wolffd@0
|
35 gdt = g*dt;
|
|
wolffd@0
|
36 ydt = y*dt;
|
|
wolffd@0
|
37 ldt = l*dt;
|
|
wolffd@0
|
38 rdt = r*dt;
|
|
wolffd@0
|
39 xdt = x*dt;
|
|
wolffd@0
|
40 [numChannels dataLength] = size(data);
|
|
wolffd@0
|
41
|
|
wolffd@0
|
42 % Initial values
|
|
wolffd@0
|
43 kt = g*A/(A+B);
|
|
wolffd@0
|
44 spont = M*y*kt/(l*kt+y*(l+r));
|
|
wolffd@0
|
45 c = spont * ones(numChannels,1);
|
|
wolffd@0
|
46 q = c*(l+r)/kt;
|
|
wolffd@0
|
47 w = c*r/x;
|
|
wolffd@0
|
48 zeroVector = zeros(numChannels,1);
|
|
wolffd@0
|
49
|
|
wolffd@0
|
50 % Now iterate through each time slice of the data. Use the
|
|
wolffd@0
|
51 % max function to implement the "if (0>" test.
|
|
wolffd@0
|
52 y = zeros(numChannels, dataLength);
|
|
wolffd@0
|
53 for i = 1:dataLength
|
|
wolffd@0
|
54 limitedSt = max(data(:,i)+A,0);
|
|
wolffd@0
|
55 kt = gdt*limitedSt./(limitedSt+B);
|
|
wolffd@0
|
56 replenish = max(ydt*(M-q),zeroVector);
|
|
wolffd@0
|
57 eject = kt.*q;
|
|
wolffd@0
|
58 loss = ldt.*c;
|
|
wolffd@0
|
59 reuptake = rdt.*c;
|
|
wolffd@0
|
60 reprocess = xdt.*w;
|
|
wolffd@0
|
61
|
|
wolffd@0
|
62 q = q + replenish - eject + reprocess;
|
|
wolffd@0
|
63 c = c + eject - loss - reuptake;
|
|
wolffd@0
|
64 w = w + reuptake - reprocess;
|
|
wolffd@0
|
65 y(:,i) = c;
|
|
wolffd@0
|
66 end
|
|
wolffd@0
|
67
|
|
wolffd@0
|
68 y = h .* y;
|
|
wolffd@0
|
69
|
|
wolffd@0
|
70 if (subtractSpont > 0)
|
|
wolffd@0
|
71 y=max(0,y-spont);
|
|
wolffd@0
|
72 end
|
