dicklyon@473: % Copyright 2012, Google, Inc.
dicklyon@473: % Author: Richard F. Lyon
dicklyon@473: %
dicklyon@473: % This Matlab file is part of an implementation of Lyon's cochlear model:
dicklyon@473: % "Cascade of Asymmetric Resonators with Fast-Acting Compression"
dicklyon@473: % to supplement Lyon's upcoming book "Human and Machine Hearing"
dicklyon@473: %
dicklyon@473: % Licensed under the Apache License, Version 2.0 (the "License");
dicklyon@473: % you may not use this file except in compliance with the License.
dicklyon@473: % You may obtain a copy of the License at
dicklyon@473: %
dicklyon@473: %     http://www.apache.org/licenses/LICENSE-2.0
dicklyon@473: %
dicklyon@473: % Unless required by applicable law or agreed to in writing, software
dicklyon@473: % distributed under the License is distributed on an "AS IS" BASIS,
dicklyon@473: % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
dicklyon@473: % See the License for the specific language governing permissions and
dicklyon@473: % limitations under the License.
dicklyon@473: 
dicklyon@473: function [zY, state] = CARFAC_CAR_Step(x_in, CAR_coeffs, state)
dicklyon@473: % function [zY, state] = CARFAC_CAR_Step(x_in, CAR_coeffs, state)
dicklyon@473: %
dicklyon@473: % One sample-time update step for the filter part of the CARFAC.
dicklyon@473: 
dicklyon@473: % Most of the update is parallel; finally we ripple inputs at the end.
dicklyon@473: 
dicklyon@473: % Local nonlinearity zA and AGC feedback zB reduce pole radius:
dicklyon@473: zA = state.zA_memory;
dicklyon@473: zB = state.zB_memory + state.dzB_memory; % AGC interpolation
dicklyon@473: r1 = CAR_coeffs.r1_coeffs;
dicklyon@473: g = state.g_memory + state.dg_memory;  % interp g
dicklyon@473: v_offset  = CAR_coeffs.v_offset;
dicklyon@473: v2_corner = CAR_coeffs.v2_corner;
dicklyon@473: v_damp_max = CAR_coeffs.v_damp_max;
dicklyon@473: 
dicklyon@473: % zB and zA are "extra damping", and multiply zr (compressed theta):
dicklyon@473: r = r1 - CAR_coeffs.zr_coeffs .* (zA + zB); 
dicklyon@473: 
dicklyon@473: % now reduce state by r and rotate with the fixed cos/sin coeffs:
dicklyon@473: z1 = r .* (CAR_coeffs.a0_coeffs .* state.z1_memory - ...
dicklyon@473:   CAR_coeffs.c0_coeffs .* state.z2_memory);
dicklyon@473: % z1 = z1 + inputs;
dicklyon@473: z2 = r .* (CAR_coeffs.c0_coeffs .* state.z1_memory + ...
dicklyon@473:   CAR_coeffs.a0_coeffs .* state.z2_memory);
dicklyon@473: 
dicklyon@473: % update the "velocity" for cubic nonlinearity, into zA:
dicklyon@473: zA = (((state.z2_memory - z2) .* CAR_coeffs.velocity_scale) + ...
dicklyon@473:   v_offset) .^ 2;
dicklyon@473: % soft saturation to make it more like an "essential" nonlinearity:
dicklyon@473: zA = v_damp_max * zA ./ (v2_corner + zA);
dicklyon@473: 
dicklyon@473: zY = CAR_coeffs.h_coeffs .* z2;  % partial output
dicklyon@473: 
dicklyon@473: % Ripple input-output path, instead of parallel, to avoid delay...
dicklyon@473: % this is the only part that doesn't get computed "in parallel":
dicklyon@473: in_out = x_in;
dicklyon@473: for ch = 1:length(zY)
dicklyon@473:   % could do this here, or later in parallel:
dicklyon@473:   z1(ch) = z1(ch) + in_out;
dicklyon@473:   % ripple, saving final channel outputs in zY
dicklyon@473:   in_out = g(ch) * (in_out + zY(ch));
dicklyon@473:   zY(ch) = in_out;
dicklyon@473: end
dicklyon@473: 
dicklyon@473: % put new state back in place of old
dicklyon@473: % (z1 and z2 are genuine temps; the others can update by reference in C)
dicklyon@473: state.z1_memory = z1;
dicklyon@473: state.z2_memory = z2;
dicklyon@473: state.zA_memory = zA;
dicklyon@473: state.zB_memory = zB;
dicklyon@473: state.zY_memory = zY;
dicklyon@473: state.g_memory = g;
dicklyon@473: