Start v3 implemention

matlab/CARFAC_Design.m

@@ -17,9 +17,10 @@
 % limitations under the License.
 
 function CF = CARFAC_Design(n_ears, fs, ...
-  CF_CAR_params, CF_AGC_params, CF_IHC_params, CF_IHC_keyword)
+  CF_CAR_params, CF_AGC_params, CF_IHC_params, CF_SYN_params, ...
+  CF_version_keyword)
 % function CF = CARFAC_Design(n_ears, fs, ...
-%   CF_CAR_params, CF_AGC_params, CF_IHC_params, CF_IHC_keyword)
+%   CF_CAR_params, CF_AGC_params, CF_IHC_params, CF_version_keyword)
 %
 % This function designs the CARFAC (Cascade of Asymmetric Resonators with
 % Fast-Acting Compression); that is, it take bundles of parameters and
@@ -32,8 +33,8 @@
 % CF_AGC_params bundles all the automatic gain control parameters
 % CF_IHC_params bundles all the inner hair cell parameters
 % Indendent of how many of these are provided, if the last arg is a string
-% it will be interpreted as an IHC_style keyword ('just_hwr' or 'one_cap';
-% omit or 'two_cap' for default v2 two-cap IHC model.
+% it will be interpreted as a CF_version_keyword ('just_hwr' or 'one_cap';
+% omit or 'two_cap' for default v2 two-cap IHC model; or 'do_syn')
 %
 % See other functions for designing and characterizing the CARFAC:
 % [naps, CF] = CARFAC_Run(CF, input_waves)
@@ -63,23 +64,30 @@
   case 5
     last_arg = CF_IHC_params;
   case 6
-    last_arg = CF_IHC_keyword;
+    last_arg = CF_SYN_params;
+  case 7
+    last_arg = CF_version_keyword;
 end
+
+
+% Last arg being a keyword can be 'do_syn' or an IHC version.
 if ischar(last_arg)  % string is a character array, not a string array.
-  IHC_keyword = last_arg;
+  CF_version_keyword = last_arg;
   num_args = nargin - 1;
 else
-  IHC_keyword = 'two_cap';  % the CARFAC v2 default
+  CF_version_keyword = 'two_cap';  % The CARFAC v2 default
   num_args = nargin;
 end
 
+% Now num_args does not count the version_keyword.  Can be up to 6.
+
 if num_args < 1
   n_ears = 1;  % if more than 1, make them identical channels;
   % then modify the design if necessary for different reasons
 end
 
 if num_args < 2
-  fs = 22050;
+  fs = 22050;  % Keeping this poor default in v3 to encourage users to always specify.
 end
 
 if num_args < 3
@@ -110,22 +118,27 @@
     'AGC_mix_coeff', 0.5);
 end
 
-if num_args < 5
-  one_cap = 0;         % bool; 1 for Allen model, 0 for new default two-cap
-  just_hwr = 0;        % bool; 0 for normal/fancy IHC; 1 for HWR
-  switch IHC_keyword
-    case 'just_hwr'
-      just_hwr = 1;        % bool; 0 for normal/fancy IHC; 1 for HWR
-    case 'one_cap'
-      one_cap = 1;         % bool; 1 for Allen model, as text states we use
-    case 'two_cap'
-      % nothing to do; accept the default
-    otherwise
-      error('unknown IHC_keyword in CARFAC_Design')
-  end
+one_cap = 0;         % bool; 1 for Allen model, 0 for new default two-cap
+just_hwr = 0;        % bool; 0 for normal/fancy IHC; 1 for HWR
+do_syn = 0;
+switch CF_version_keyword
+  case 'just_hwr'
+    just_hwr = 1;        % bool; 0 for normal/fancy IHC; 1 for HWR
+  case 'one_cap'
+    one_cap = 1;         % bool; 1 for Allen model, as text states we use
+  case 'do_syn';
+    do_syn = 1;
+  case 'two_cap'
+    % nothing to do; accept the v2 default, two-cap IHC, no SYN.
+  otherwise
+    error('unknown IHC_keyword in CARFAC_Design')
+end
+
+if num_args < 5  % Default IHC_params
   CF_IHC_params = struct( ...
     'just_hwr', just_hwr, ...  % not just a simple HWR
     'one_cap', one_cap, ...    % bool; 0 for new two-cap hack
+    'do_syn', do_syn, ...      % bool; 1 for v3 synapse feature
     'tau_lpf', 0.000080, ...   % 80 microseconds smoothing twice
     'tau_out', 0.0005, ...     % depletion tau is pretty fast
     'tau_in', 0.010, ...       % recovery tau is slower
@@ -135,6 +148,26 @@
     'tau2_in', 0.010);         % recovery tau is slower 10 ms
 end
 
+if num_args < 6  % Default the SYN_params.
+  n_classes = 3;  % Default.  Modify params and redesign to change.
+  IHCs_per_channel = 10;  % Maybe 20 would be better?
+  % Parameters could generally have columns if class-dependent.
+  CF_SYN_params = struct( ...
+    'do_syn', do_syn, ...  % This may just turn it off completely.
+    'fs', fs, ...
+    'n_classes', n_classes, ...
+    'healthy_n_fibers', [50, 35, 25] * IHCs_per_channel, ... 
+    'spont_rates', [50, 6, 1], ...
+    'sat_rates', 200, ...
+    'sat_reservoir', 0.2, ...
+    'v_width', 0.02, ...
+    'tau_lpf', 0.000080, ...
+    'reservoir_tau', 0.020, ...
+    'agc_weights', [1.2, 1.2, 1.2] / (22050*IHCs_per_channel));  % Tweaked.
+  % The weights 1.2 were picked before correctly account for sample rate
+  % and number of fibers.  This way works for more different numbers.
+end
+
 % first figure out how many filter stages (PZFC/CARFAC channels):
 pole_Hz = CF_CAR_params.first_pole_theta * fs / (2*pi);
 n_ch = 0;
@@ -160,13 +193,15 @@
 CAR_coeffs = CARFAC_DesignFilters(CF_CAR_params, fs, pole_freqs);
 AGC_coeffs = CARFAC_DesignAGC(CF_AGC_params, fs, n_ch);
 IHC_coeffs = CARFAC_DesignIHC(CF_IHC_params, fs, n_ch);
+SYN_coeffs = CARFAC_DesignSynapses(CF_SYN_params, fs, pole_freqs);
 
 % Copy same designed coeffs into each ear (can do differently in the
 % future, e.g. for unmatched OHC_health).
 for ear = 1:n_ears
   ears(ear).CAR_coeffs = CAR_coeffs;
   ears(ear).AGC_coeffs = AGC_coeffs;
   ears(ear).IHC_coeffs = IHC_coeffs;
+  ears(ear).SYN_coeffs = SYN_coeffs;
 end
 
 CF = struct( ...
@@ -175,12 +210,14 @@
   'CAR_params', CF_CAR_params, ...
   'AGC_params', CF_AGC_params, ...
   'IHC_params', CF_IHC_params, ...
+  'SYN_params', CF_SYN_params, ...
   'n_ch', n_ch, ...
   'pole_freqs', pole_freqs, ...
   'ears', ears, ...
   'n_ears', n_ears, ...
   'open_loop', 0, ...
-  'linear_car', 0);
+  'linear_car', 0, ...
+  'do_syn', do_syn);
 
 
 %% Design the filter coeffs:
@@ -539,3 +576,67 @@
       'ihc_accum', 0);
   end
 end
+
+%% the SYN design coeffs:
+function SYN_coeffs = CARFAC_DesignSynapses(SYN_params, fs, pole_freqs)
+
+if ~SYN_params.do_syn
+  SYN_coeffs = [];  % Just return empty if we're not doing SYN.
+  return
+end
+
+n_ch = length(pole_freqs);
+n_classes = SYN_params.n_classes;
+
+v_widths = SYN_params.v_width * ones(1, n_classes);
+
+% Do some design.  First, gains to get sat_rate when sigmoid is 1, which
+% involves the reservoir steady-state solution.
+% Most of these are not per-channel, but could be expanded that way
+% later if desired.
+
+% Mean sat prob of spike per sample per neuron, likely same for all
+% classes.
+% Use names 1 for sat and 0 for spont in some of these.
+p1 = SYN_params.sat_rates / fs;
+p0 = SYN_params.spont_rates / fs;
+
+w1 = SYN_params.sat_reservoir;
+q1 = 1 - w1;
+% Assume the sigmoid is switching between 0 and 1 at 50% duty cycle, so
+% normalized mean value is 0.5 at saturation.
+s1 = 0.5;
+r1 = s1*w1;
+% solve q1 = a1*r1 for gain coeff a1:
+a1 = q1 ./ r1;
+% solve p1 = a2*r1 for gain coeff a2:
+a2 = p1 ./ r1;
+
+% Now work out how to get the desired spont.
+r0 = p0 ./ a2;
+q0 = r0 .* a1;
+w0 = 1 - q0;
+s0 = r0 ./ w0;
+% Solve for (negative) sigmoid midpoint offset that gets s0 right.
+offsets = log((1 - s0)./s0);
+
+spont_p = a2 .* w0 .* s0;  % should match p0; check it; yes it does.
+
+agc_weights = fs * SYN_params.agc_weights;
+spont_sub = (SYN_params.healthy_n_fibers .* spont_p) * agc_weights';
+
+% Copy stuff needed at run time into coeffs.
+SYN_coeffs = struct( ...
+  'n_ch', n_ch, ...
+  'n_classes', n_classes, ...
+  'n_fibers', ones(n_ch,1) * SYN_params.healthy_n_fibers, ...
+  'v_widths', v_widths, ...
+  'v_halfs', offsets .* v_widths, ...  % Same units as v_recep and v_widths.
+  'a1', a1, ...  % Feedback gain
+  'a2', a2, ...  % Output gain
+  'agc_weights', agc_weights, ... % For making a nap out to agc in.
+  'spont_p', spont_p, ... % used only to init the output LPF
+  'spont_sub', spont_sub, ...
+  'res_lpf_inits', q0, ...
+  'res_coeff', 1 - exp(-1/(SYN_params.reservoir_tau * fs)), ...
+  'lpf_coeff', 1 - exp(-1/(SYN_params.tau_lpf * fs)));

matlab/CARFAC_IHC_Step.m

@@ -16,12 +16,16 @@
 % See the License for the specific language governing permissions and
 % limitations under the License.
 
-function [ihc_out, state] = CARFAC_IHC_Step(bm_out, coeffs, state);
+function [ihc_out, state, v_recep] = CARFAC_IHC_Step(bm_out, coeffs, state);
 % function [ihc_out, state] = CARFAC_IHC_Step(bm_out, coeffs, state);
 %
 % One sample-time update of inner-hair-cell (IHC) model, including the
 % detection nonlinearity and one or two capacitor state variables.
+%
+% receptor_potential output will be empty except in two_cap mode.  
+% Use it as input to CARFAC_SYN_Step to model synapses to get firing rates.
 
+v_recep = [];  % For cases other than two_cap and do_syn it's not used.
 if coeffs.just_hwr
   ihc_out = min(2, max(0, bm_out));  % limit it for stability
 else
@@ -63,7 +67,10 @@
     state.lpf1_state = state.lpf1_state + coeffs.lpf_coeff * ...
       (ihc_out - state.lpf1_state);
     ihc_out = state.lpf1_state - coeffs.rest_output;
+    % Return a modified receptor potential that's zero at rest, for SYN.
+    v_recep = coeffs.rest_cap1 - state.cap1_voltage;
   end
 end
 
+% Leaving this for v2 cochleagram compatibility, but no v3/SYN version:
 state.ihc_accum = state.ihc_accum + ihc_out;  % for where decimated output is useful

matlab/CARFAC_Init.m

@@ -29,6 +29,9 @@
   CF.ears(ear).CAR_state = CAR_Init_State(CF.ears(ear).CAR_coeffs);
   CF.ears(ear).IHC_state = IHC_Init_State(CF.ears(ear).IHC_coeffs);
   CF.ears(ear).AGC_state = AGC_Init_State(CF.ears(ear).AGC_coeffs);
+  if CF.do_syn
+     CF.ears(ear).SYN_state = SYN_Init_State(CF.ears(ear).SYN_coeffs);
+  end
 end
 
 
@@ -80,3 +83,11 @@
       );
   end
 end
+
+function state = SYN_Init_State(coeffs)
+n_ch = coeffs.n_ch;
+n_cl = coeffs.n_classes;
+state = struct( ...
+  'reservoirs', ones(n_ch, 1) * coeffs.res_lpf_inits, ...  % 0 full, 1 empty.
+  'lpf_state', ones(n_ch, 1) * coeffs.spont_p); 
+

matlab/CARFAC_Run_Segment.m

@@ -111,15 +111,30 @@
       input_waves(k, ear), CF.ears(ear).CAR_coeffs, CF.ears(ear).CAR_state);
 
     % update IHC state & output on every time step, too
-    [ihc_out, CF.ears(ear).IHC_state] = CARFAC_IHC_Step( ...
+    [ihc_out, CF.ears(ear).IHC_state, v_recep] = CARFAC_IHC_Step( ...
       car_out, CF.ears(ear).IHC_coeffs, CF.ears(ear).IHC_state);
 
-    % run the AGC update step, decimating internally,
+    if CF.do_syn
+      % Use v_recep from IHC_Step to
+      % update the SYNapse state and get firings and new nap.
+      [syn_out, firings, CF.ears(ear).SYN_state] = CARFAC_SYN_Step( ...
+        v_recep, CF.ears(ear).SYN_coeffs, CF.ears(ear).SYN_state);
+      % Use sum over syn_outs classes, appropriately scaled, as nap to agc.
+      % firings always positive, unless v2 ihc_out.
+      % syn_out can go a little negative; should be zero at rest.
+      nap = syn_out;
+      % Maybe still should add a way to return firings (of the classes).
+    else
+      % v2, ihc_out already has rest_output subtracted.
+      nap = ihc_out;  % If no SYN, ihc_out goes to nap and agc as in v2.
+    end
+
+    % Use nap to run the AGC update step, decimating internally,
     [CF.ears(ear).AGC_state, AGC_updated] = CARFAC_AGC_Step( ...
-      ihc_out, CF.ears(ear).AGC_coeffs, CF.ears(ear).AGC_state);
+      nap, CF.ears(ear).AGC_coeffs, CF.ears(ear).AGC_state);
 
     % save some output data:
-    naps(k, :, ear) = ihc_out;  % output to neural activity pattern
+    naps(k, :, ear) = nap;  % output to neural activity pattern
     if do_BM
       BM(k, :, ear) = car_out;
       state = CF.ears(ear).CAR_state;

matlab/CARFAC_SYN_Step.m

@@ -0,0 +1,31 @@
+function [syn_out, firings, state] = CARFAC_SYN_Step(v_recep, coeffs, state)
+% Drive multiple synapse classes with receptor potential from IHC,
+% returning instantaneous spike rates per class, for a group of neurons
+% associated with the CF channel, including reductions due to synaptopathy.
+
+% Normalized offset position into neurotransmitter release sigmoid.
+x = (v_recep - coeffs.v_halfs) ./ coeffs.v_widths;
+
+s = 1 ./ (1 + exp(-x));  % Between 0 and 1; positive at rest.
+q = state.reservoirs;  % aka 1 - w, between 0 and 1; positive at rest.
+r = (1 - q) .* s;  % aka w*s, between 0 and 1, proportional to release rate.
+
+% Smooth once with LPF (receptor potential was already smooth), after
+% applying the gain coeff a2 to convert to firing prob per sample.
+state.lpf_state = state.lpf_state + coeffs.lpf_coeff * ...
+  (coeffs.a2 .* r - state.lpf_state);  % this is firing probs.
+firing_probs = state.lpf_state;  % Poisson rate per neuron per sample.
+% Include number of effective neurons per channel here, and interval T;
+% so the rates (instantaneous action potentials per second) can be huge.
+firings = coeffs.n_fibers .* firing_probs;
+
+% Feedback, to update reservoir state q for next time.
+state.reservoirs = q + coeffs.res_coeff .* (coeffs.a1 .* r - q);
+
+% Make an output that resembles ihc_out, to go to agc_in (collapse over classes).
+% Includes synaptopathy's presumed effect of reducing feedback via n_fibers.
+% But it's relative to the healthy nominal spont, so could potentially go
+% a bit negative in quiet is there was loss of high-spont or medium-spont units.
+% The weight multiplication is an inner product, reducing n_classes
+% columns to 1 column (first transpose the agc_weights row to a column).
+syn_out = (coeffs.n_fibers .* firing_probs) * coeffs.agc_weights' - coeffs.spont_sub;