Start v3 implemention

CARFAC_IHC_Step.m

@@ -0,0 +1,75 @@
+% Copyright 2012 The CARFAC Authors. All Rights Reserved.
+% Author: Richard F. Lyon
+%
+% This file is part of an implementation of Lyon's cochlear model:
+% "Cascade of Asymmetric Resonators with Fast-Acting Compression"
+%
+% Licensed under the Apache License, Version 2.0 (the "License");
+% you may not use this file except in compliance with the License.
+% You may obtain a copy of the License at
+%
+%     http://www.apache.org/licenses/LICENSE-2.0
+%
+% 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 [ihc_out, state, receptor_potential] = 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.
+
+if coeffs.just_hwr
+  ihc_out = min(2, max(0, bm_out));  % limit it for stability
+  receptor_potential = [];
+else
+  conductance = CARFAC_Detect(bm_out);  % rectifying nonlinearity
+
+  if coeffs.one_cap;
+    % Output comes from receptor current like in Hall and Allen's models.
+    ihc_out = conductance .* state.cap_voltage;
+    state.cap_voltage = state.cap_voltage - ...
+      ihc_out .* coeffs.out_rate + ...
+      (1 - state.cap_voltage) .* coeffs.in_rate;
+    ihc_out = ihc_out * coeffs.output_gain;
+    % Smooth it twice with LPF:
+    state.lpf1_state = state.lpf1_state + coeffs.lpf_coeff * ...
+      (ihc_out - state.lpf1_state);
+    state.lpf2_state = state.lpf2_state + coeffs.lpf_coeff * ...
+      (state.lpf1_state - state.lpf2_state);
+    ihc_out = state.lpf2_state - coeffs.rest_output;
+    receptor_potential = [];
+  else
+    % Change to 2-cap version mediated by receptor potential at cap1:
+    % Geisler book fig 8.4 suggests 40 to 800 Hz corner.
+    receptor_current = conductance .* state.cap1_voltage;
+    % "out" means charge depletion; "in" means restoration toward 1.
+    state.cap1_voltage = state.cap1_voltage - ...
+      receptor_current .* coeffs.out1_rate + ...
+      (1 - state.cap1_voltage) .* coeffs.in1_rate;
+    % Amount of depletion below 1 is receptor potential.
+    receptor_potential = 1 - state.cap1_voltage;  % Already smooth.
+    % Identity map from receptor potential to neurotransmitter conductance.
+    ihc_out = receptor_potential .* state.cap2_voltage;  % Now a current.
+    % cap2 represents transmitter store; another adaptive gain.
+    % Deplete the transmitter store like in Meddis models:
+    state.cap2_voltage = state.cap2_voltage - ...
+      ihc_out .* coeffs.out2_rate + ...
+      (1 - state.cap2_voltage) .* coeffs.in2_rate;
+    % Awkwardly, gain needs to be here for the init states to be right.
+    ihc_out = ihc_out * coeffs.output_gain;
+    % smooth once more with LPF (receptor potential was already smooth):
+    state.lpf1_state = state.lpf1_state + coeffs.lpf_coeff * ...
+      (ihc_out - state.lpf1_state);
+    ihc_out = state.lpf1_state - coeffs.rest_output;
+  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

CARFAC_Init.m

@@ -0,0 +1,94 @@
+% Copyright 2012 The CARFAC Authors. All Rights Reserved.
+% Author: Richard F. Lyon
+%
+% This file is part of an implementation of Lyon's cochlear model:
+% "Cascade of Asymmetric Resonators with Fast-Acting Compression"
+%
+% Licensed under the Apache License, Version 2.0 (the "License");
+% you may not use this file except in compliance with the License.
+% You may obtain a copy of the License at
+%
+%     http://www.apache.org/licenses/LICENSE-2.0
+%
+% 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_Init(CF)
+% function CF = CARFAC_Init(CF)
+%
+% Initialize state for one or more ears of CF.
+% This allocates and zeros all the state vector storage in the CF struct.
+
+n_ears = CF.n_ears;
+
+for ear = 1:n_ears
+  % for now there's only one coeffs, not one per ear
+  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
+
+
+function state = CAR_Init_State(coeffs)
+n_ch = coeffs.n_ch;
+state = struct( ...
+  'z1_memory', zeros(n_ch, 1), ...
+  'z2_memory', zeros(n_ch, 1), ...
+  'zA_memory', zeros(n_ch, 1), ...
+  'zB_memory', coeffs.zr_coeffs, ...
+  'dzB_memory', zeros(n_ch, 1), ...
+  'zY_memory', zeros(n_ch, 1), ...
+  'g_memory', coeffs.g0_coeffs, ...
+  'dg_memory', zeros(n_ch, 1), ...
+  'ac_coupler', zeros(n_ch, 1) ...
+  );
+
+
+function state = AGC_Init_State(coeffs)
+n_ch = coeffs(1).n_ch;
+n_AGC_stages = coeffs(1).n_AGC_stages;
+state = struct([]);
+for stage = 1:n_AGC_stages
+  % Initialize state recursively...
+  state(stage).AGC_memory = zeros(n_ch, 1);
+  state(stage).input_accum = zeros(n_ch, 1);
+  state(stage).decim_phase = 0;  % integer decimator phase
+end
+
+
+function state = IHC_Init_State(coeffs)
+n_ch = coeffs.n_ch;
+if coeffs.just_hwr
+  state = struct('ihc_accum', zeros(n_ch, 1));
+else
+  if coeffs.one_cap
+    state = struct( ...
+      'ihc_accum', zeros(n_ch, 1), ...
+      'cap_voltage', coeffs.rest_cap * ones(n_ch, 1), ...
+      'lpf1_state', coeffs.rest_output * ones(n_ch, 1), ...
+      'lpf2_state', coeffs.rest_output * ones(n_ch, 1) ...
+      );
+  else
+    state = struct( ...
+      'ihc_accum', zeros(n_ch, 1), ...
+      'cap1_voltage', coeffs.rest_cap1 * ones(n_ch, 1), ...
+      'cap2_voltage', coeffs.rest_cap2 * ones(n_ch, 1), ...
+      'lpf1_state', coeffs.rest_output * ones(n_ch, 1) ...
+      );
+  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_inits, ...  % 1 means full.
+  'lpf_state', zeros(n_ch, n_cl) ); 
+% Need to make different rest fullnesses for the different classes!!!
+

CARFAC_Run_Segment.m

@@ -0,0 +1,157 @@
+% Copyright 2012 The CARFAC Authors. All Rights Reserved.
+% Author Richard F. Lyon
+%
+% This file is part of an implementation of Lyon's cochlear model:
+% "Cascade of Asymmetric Resonators with Fast-Acting Compression"
+%
+% Licensed under the Apache License, Version 2.0 (the "License");
+% you may not use this file except in compliance with the License.
+% You may obtain a copy of the License at
+%
+%     http://www.apache.org/licenses/LICENSE-2.0
+%
+% 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 [naps, CF, BM, seg_ohc, seg_agc] = CARFAC_Run_Segment(...
+  CF, input_waves)
+% function [naps, CF, BM, seg_ohc, seg_agc] = CARFAC_Run_Segment(...
+%   CF, input_waves)
+%
+% This function runs the CARFAC; that is, filters a 1 or more channel
+% sound input segment to make one or more neural activity patterns (naps);
+% it can be called multiple times for successive segments of any length,
+% as long as the returned CF with modified state is passed back in each
+% time.
+%
+% input_waves is a column vector if there's just one audio channel;
+% more generally, it has a row per time sample, a column per audio channel.
+%
+% naps has a row per time sample, a column per filterbank channel, and
+% a layer per audio channel if more than 1.
+% BM is basilar membrane motion (filter outputs before detection).
+%
+% the input_waves are assumed to be sampled at the same rate as the
+% CARFAC is designed for; a resampling may be needed before calling this.
+%
+% The function works as an outer iteration on time, updating all the
+% filters and AGC states concurrently, so that the different channels can
+% interact easily.  The inner loops are over filterbank channels, and
+% this level should be kept efficient.
+%
+% seg_ohc seg_agc are optional extra outputs useful for seeing what the
+% ohc nonlinearity and agc are doing; both in terms of extra damping.
+
+if nargout > 2
+  do_BM = 1;
+else
+  do_BM = 0;
+end
+
+[n_samp, n_ears] = size(input_waves);
+
+if n_ears ~= CF.n_ears
+  error('bad number of input_waves channels passed to CARFAC_Run')
+end
+
+if ~isfield(CF, 'open_loop')  % Find open_loop in CF or default it.
+  CF.open_loop = 0;
+end
+
+if ~isfield(CF, 'linear_car')  % Find linear in CF or default it.
+  CF.linear_car = 0;
+end
+
+if ~isfield(CF, 'use_delay_buffer')  % To let CAR be fully parallel.
+  CF.use_delay_buffer = 0;
+end
+
+n_ch = CF.n_ch;
+naps = zeros(n_samp, n_ch, n_ears);  % allocate space for result
+if do_BM
+  BM = zeros(n_samp, n_ch, n_ears);
+  seg_ohc = zeros(n_samp, n_ch, n_ears);
+  seg_agc = zeros(n_samp, n_ch, n_ears);
+end
+
+% A 2022 addition to make open-loop running behave.  In open_loop mode, these
+% coefficients are set, per AGC filter outputs, when we CARFAC_Close_AGC_Loop
+% on AGC filter output updates.  They drive the zB and g coefficients to the
+% intended value by the next update, but if open_loop they would just keep
+% going, extrapolating.  The point of open_loop mode is to stop them moving,
+% so we need to make sure these deltas are zeroed in case the mode is switched
+% from closed to open, as in some tests to evaluate the transfer functions
+% before and after adapting to a signal.
+if CF.open_loop
+  % The interpolators may be running if it was previously run closed loop.
+  for ear = 1:CF.n_ears
+    CF.ears(ear).CAR_state.dzB_memory = 0;  % To stop intepolating zB.
+    CF.ears(ear).CAR_state.dg_memory = 0;  % To stop intepolating g.
+  end
+end
+
+% Apply control coeffs to where they are needed.
+for ear = 1:CF.n_ears
+  CF.ears(ear).CAR_coeffs.linear = CF.linear_car;  % Skip OHC nonlinearity.
+  CF.ears(ear).CAR_coeffs.use_delay_buffer = CF.use_delay_buffer;
+end
+
+detects = zeros(n_ch, n_ears);
+for k = 1:n_samp
+  % at each time step, possibly handle multiple channels
+  for ear = 1:n_ears
+
+    % This would be cleaner if we could just get and use a reference to
+    % CF.ears(ear), but Matlab doesn't work that way...
+
+    [car_out, CF.ears(ear).CAR_state] = CARFAC_CAR_Step( ...
+      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, v_recep] = CARFAC_IHC_Step( ...
+      car_out, CF.ears(ear).IHC_coeffs, CF.ears(ear).IHC_state);
+
+    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( ...
+      nap, CF.ears(ear).AGC_coeffs, CF.ears(ear).AGC_state);
+
+    % save some output data:
+    naps(k, :, ear) = nap;  % output to neural activity pattern
+    if do_BM
+      BM(k, :, ear) = car_out;
+      state = CF.ears(ear).CAR_state;
+      seg_ohc(k, :, ear) = state.zA_memory;
+      seg_agc(k, :, ear) = state.zB_memory;
+    end
+  end
+
+  % connect the feedback from AGC_state to CAR_state when it updates;
+  % all ears together here due to mixing across them:
+  if AGC_updated
+    if n_ears > 1
+      % do multi-aural cross-coupling:
+      CF.ears = CARFAC_Cross_Couple(CF.ears);
+    end
+    if ~CF.open_loop
+      CF = CARFAC_Close_AGC_Loop(CF);  % Starts the interpolation of zB and g.
+    end
+  end
+end

CARFAC_SYN_Step.m

@@ -0,0 +1,57 @@
+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.
+
+if nargin < 1
+  error('SYN_step needs at least a v_recep input.')
+end
+
+if nargin < 3  % Temporarily do the setup defaulting here; design it elsewhere later.
+  n_ones = ones(1, 3);
+  v_width = 0.05;
+  % Parameters could generally have columns if channel-dependent.
+  max_rate = 2500;  % % Instantaneous max at onset, per Kiang figure 5.18.
+  fs = 48000;
+  coeffs = struct( ...
+    'fs', fs, ...
+    'n_classes', length(n_ones), ...
+    'max_probs', (max_rate / fs) * n_ones, ...  
+    'n_fibers', [50, 35, 25], ...  % Synaptopathy comes in here; channelize it, too.
+    'v_width', v_width * n_ones, ...  % Try to match IHC out units.
+    'v_half', v_width .* [3, 5, 7], ...  % Same units as v_width and v_recep.
+    'out_rate', 0.02, ...  % Depletion can be quick (few ms).
+    'recovery', 1e-4);  % Recovery tau about 10,000 sample times.
+end
+
+if nargin < 2  % Init the state, sized according to the coeffs.
+  n_ch = 1;
+  % "Full" reservoir state is 1.0, empty 0.0.
+  state = struct( ...
+    'reservoirs', ones(n_ch, coeffs.n_classes));
+end
+
+% This release_rates is relative to a max of 1, usually way lower.
+release_rates = state.reservoirs ./ ...
+  (1 + exp(-(v_recep - coeffs.v_half)./coeffs.v_width));
+
+% Smooth once with LPF (receptor potential was already smooth):
+state.lpf_state = state.lpf_state + coeffs.lpf_coeff * ...
+  (release_rates - state.lpf_state);
+release_rates = state.lpf_state;
+
+% Deplete reservoirs some, independent of synaptopathy numbers (we assume).
+state.reservoirs = state.reservoirs - coeffs.out_rate .* release_rates;
+% Refill reservoirs a little toward 1.
+state.reservoirs = state.reservoirs + coeffs.recovery .* (1 - state.reservoirs);
+
+% Includes 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 .* coeffs.max_probs .* release_rates;
+
+% Make an output that resembles ihc_out, to go the agc_in.
+syn_out = firings * coeffs.agc_weights_col - coeffs.rest_output;
+
+
+
+