Copybara import of the project:

--
628c358 by Richard Lyon <dicklyon@google.com>:

Start v3 implemention

Adding new SYN stage (synape model) for a v3 CARFAC.
--
3ab484e by Richard Lyon <dicklyon@google.com>:

Add files via upload
--
ab450a0 by Richard Lyon <dicklyon@google.com>:

Fixes per robsc's comments on PR 18

Fixed Design process per robsc's comments on PR 18.

--
6975bde by Richard Lyon <dicklyon@google.com>:

Some SYN improvements; still needs a lot of work.

Some SYN improvements, but really it needs a revamp of the Design process, with better parameterization and automatic calculation of rest rate to subtract to get AGC back closer to how it worked in old versions.

--
17f5564 by Richard Lyon <dicklyon@google.com>:

Reparameterized, better design phase

I reparameterized the SYN stage in terms of spontaneous and saturation rates for the classes and few other things, and implemented the design process to account for everything to get those to come out right, all in terms of variables that show up on the diagram that's an augmentation of the v2 two_cap IHC diagram (for documentation, in progress).
Still need to move some of the numbers into SYN_params.

--
b79778a by Richard Lyon <dicklyon@google.com>:

Move SYN constants from code into default parameters

Moving a bunch of magic numbers out of the SYN_Design code and into the default SYN_params, so they are modifiable by user code.  And various other code tweaks and naming improvements (e.g. changed z to r for release rate; short variable names to facilitate labeling on doc diagrams).

--
0181d1e by Richard Lyon <dicklyon@google.com>:

Make agc_weights work better across sample rate and n_fiber changes

The agc_weights needed to be appropriately set up to do the right thing with different sample rates, and to allow different numbers of IHCs per channel.  The AGC feedback is still proportional to the number of fibers, healthy or not, so the agc_weights parameters might need adjustment if the healthy_n_fibers is changed.  The param n_fibers was renamed healthy_n_fibers, as that's what applied at design time, and as the n_fibers in the coeffs may change with model training (in the JAX version at least).

--
6cbf240 by Richard Lyon <dicklyon@google.com>:

Delete .DS_Store

COPYBARA_INTEGRATE_REVIEW=#18 from google:carfac-v3-start 6cbf240
PiperOrigin-RevId: 681626192

matlab/CARFAC_Design.m

@@ -9,17 +9,18 @@
 % You may obtain a copy of the License at
 %
 %     http://www.apache.org/licenses/LICENSE-2.0
-%
+%l
 % Unless required by applicable law or agreed to in writing, software
 % distributed under the License is distributed on an "AS IS" BASIS,
 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 % See the License for the specific language governing permissions and
 % 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;