alexbrandmeyer@620: //
alexbrandmeyer@620: //  ear.cc
alexbrandmeyer@620: //  CARFAC Open Source C++ Library
alexbrandmeyer@620: //
alexbrandmeyer@620: //  Created by Alex Brandmeyer on 5/10/13.
alexbrandmeyer@620: //
alexbrandmeyer@620: // This C++ file is part of an implementation of Lyon's cochlear model:
alexbrandmeyer@620: // "Cascade of Asymmetric Resonators with Fast-Acting Compression"
alexbrandmeyer@620: // to supplement Lyon's upcoming book "Human and Machine Hearing"
alexbrandmeyer@620: //
alexbrandmeyer@620: // Licensed under the Apache License, Version 2.0 (the "License");
alexbrandmeyer@620: // you may not use this file except in compliance with the License.
alexbrandmeyer@620: // You may obtain a copy of the License at
alexbrandmeyer@620: //
alexbrandmeyer@620: //     http://www.apache.org/licenses/LICENSE-2.0
alexbrandmeyer@620: //
alexbrandmeyer@620: // Unless required by applicable law or agreed to in writing, software
alexbrandmeyer@620: // distributed under the License is distributed on an "AS IS" BASIS,
alexbrandmeyer@620: // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
alexbrandmeyer@620: // See the License for the specific language governing permissions and
alexbrandmeyer@620: // limitations under the License.
alexbrandmeyer@620: 
ronw@670: #include <assert.h>
ronw@670: 
alexbrandmeyer@620: #include "ear.h"
alexbrandmeyer@620: 
alexbrandmeyer@621: // The 'InitEar' function takes a set of model parameters and initializes the
alexbrandmeyer@621: // design coefficients and model state variables needed for running the model
alexbrandmeyer@668: // on a single audio channel.
alexbrandmeyer@668: void Ear::InitEar(const int n_ch, const FPType fs,
alexbrandmeyer@668:                   const FloatArray& pole_freqs, const CARParams& car_params,
alexbrandmeyer@668:                   const IHCParams& ihc_params, const AGCParams& agc_params) {
alexbrandmeyer@621:   // The first section of code determines the number of channels that will be
alexbrandmeyer@621:   // used in the model on the basis of the sample rate and the CAR parameters
alexbrandmeyer@621:   // that have been passed to this function.
alexbrandmeyer@621:   n_ch_ = n_ch;
alexbrandmeyer@621:   // These functions use the parameters that have been passed to design the
alexbrandmeyer@668:   // coefficients for the first two stages of the model.
alexbrandmeyer@668:   car_coeffs_.Design(car_params, fs, pole_freqs);
alexbrandmeyer@668:   ihc_coeffs_.Design(ihc_params, fs);
alexbrandmeyer@668:   // This code initializes the coefficients for each of the AGC stages.
alexbrandmeyer@668:   agc_coeffs_.resize(agc_params.n_stages_);
alexbrandmeyer@668:   FPType previous_stage_gain = 0.0;
alexbrandmeyer@668:   FPType decim = 1.0;
alexbrandmeyer@668:   for (int stage = 0; stage < agc_params.n_stages_; ++stage) {
alexbrandmeyer@668:     agc_coeffs_[stage].Design(agc_params, stage, fs, previous_stage_gain,
alexbrandmeyer@668:                               decim);
alexbrandmeyer@668:     // Two variables store the decimation and gain levels for use in the design
alexbrandmeyer@668:     // of the subsequent stages.
alexbrandmeyer@668:     previous_stage_gain = agc_coeffs_[stage].agc_gain_;
alexbrandmeyer@668:     decim = agc_coeffs_[stage].decim_;
alexbrandmeyer@668:   }
alexbrandmeyer@621:   // Once the coefficients have been determined, we can initialize the state
alexbrandmeyer@621:   // variables that will be used during runtime.
alexbrandmeyer@621:   InitCARState();
alexbrandmeyer@621:   InitIHCState();
alexbrandmeyer@621:   InitAGCState();
alexbrandmeyer@620: }
alexbrandmeyer@620: 
alexbrandmeyer@621: void Ear::InitCARState() {
alexbrandmeyer@621:   car_state_.z1_memory_ = FloatArray::Zero(n_ch_);
alexbrandmeyer@621:   car_state_.z2_memory_ = FloatArray::Zero(n_ch_);
alexbrandmeyer@621:   car_state_.za_memory_ = FloatArray::Zero(n_ch_);
alexbrandmeyer@621:   car_state_.zb_memory_ = car_coeffs_.zr_coeffs_;
alexbrandmeyer@621:   car_state_.dzb_memory_ = FloatArray::Zero(n_ch_);
alexbrandmeyer@621:   car_state_.zy_memory_ = FloatArray::Zero(n_ch_);
alexbrandmeyer@621:   car_state_.g_memory_ = car_coeffs_.g0_coeffs_;
alexbrandmeyer@621:   car_state_.dg_memory_ = FloatArray::Zero(n_ch_);
alexbrandmeyer@621: }
alexbrandmeyer@621: 
alexbrandmeyer@621: void Ear::InitIHCState() {
alexbrandmeyer@621:   ihc_state_.ihc_accum_ = FloatArray::Zero(n_ch_);
alexbrandmeyer@621:   if (! ihc_coeffs_.just_hwr_) {
alexbrandmeyer@621:     ihc_state_.ac_coupler_ = FloatArray::Zero(n_ch_);
alexbrandmeyer@621:     ihc_state_.lpf1_state_ = ihc_coeffs_.rest_output_ * FloatArray::Ones(n_ch_);
alexbrandmeyer@621:     ihc_state_.lpf2_state_ = ihc_coeffs_.rest_output_ * FloatArray::Ones(n_ch_);
alexbrandmeyer@621:     if (ihc_coeffs_.one_cap_) {
alexbrandmeyer@621:       ihc_state_.cap1_voltage_ = ihc_coeffs_.rest_cap1_ *
alexbrandmeyer@668:       FloatArray::Ones(n_ch_);
alexbrandmeyer@621:     } else {
alexbrandmeyer@621:       ihc_state_.cap1_voltage_ = ihc_coeffs_.rest_cap1_ *
alexbrandmeyer@668:       FloatArray::Ones(n_ch_);
alexbrandmeyer@621:       ihc_state_.cap2_voltage_ = ihc_coeffs_.rest_cap2_ *
alexbrandmeyer@668:       FloatArray::Ones(n_ch_);
alexbrandmeyer@621:     }
alexbrandmeyer@621:   }
alexbrandmeyer@621: }
alexbrandmeyer@621: 
alexbrandmeyer@621: void Ear::InitAGCState() {
alexbrandmeyer@668:   int n_agc_stages = agc_coeffs_.size();
alexbrandmeyer@668:   agc_state_.resize(n_agc_stages);
alexbrandmeyer@668:   for (int i = 0; i < n_agc_stages; ++i) {
alexbrandmeyer@668:     agc_state_[i].decim_phase_ = 0;
alexbrandmeyer@668:     agc_state_[i].agc_memory_ = FloatArray::Zero(n_ch_);
alexbrandmeyer@668:     agc_state_[i].input_accum_ = FloatArray::Zero(n_ch_);
alexbrandmeyer@621:   }
alexbrandmeyer@621: }
alexbrandmeyer@621: 
alexbrandmeyer@668: void Ear::CARStep(const FPType input, FloatArray* car_out) {
alexbrandmeyer@668:   // This interpolates g.
alexbrandmeyer@668:   car_state_.g_memory_ = car_state_.g_memory_ + car_state_.dg_memory_;
alexbrandmeyer@621:   // This calculates the AGC interpolation state.
alexbrandmeyer@668:   car_state_.zb_memory_ = car_state_.zb_memory_ + car_state_.dzb_memory_;
alexbrandmeyer@621:   // This updates the nonlinear function of 'velocity' along with zA, which is
alexbrandmeyer@621:   // a delay of z2.
alexbrandmeyer@668:   FloatArray nonlinear_fun(n_ch_);
alexbrandmeyer@668:   FloatArray velocities = car_state_.z2_memory_ - car_state_.za_memory_;
alexbrandmeyer@668:   OHCNonlinearFunction(velocities, &nonlinear_fun);
alexbrandmeyer@668:   // Here, zb_memory_ * nonlinear_fun is "undamping" delta r.
alexbrandmeyer@668:   FloatArray r = car_coeffs_.r1_coeffs_ + (car_state_.zb_memory_ *
alexbrandmeyer@668:                                            nonlinear_fun);
alexbrandmeyer@668:   car_state_.za_memory_ = car_state_.z2_memory_;
alexbrandmeyer@621:   // Here we reduce the CAR state by r and rotate with the fixed cos/sin coeffs.
alexbrandmeyer@668:   FloatArray z1  = r * ((car_coeffs_.a0_coeffs_ * car_state_.z1_memory_) -
alexbrandmeyer@668:                         (car_coeffs_.c0_coeffs_ * car_state_.z2_memory_));
alexbrandmeyer@668:   car_state_.z2_memory_ = r *
alexbrandmeyer@668:   ((car_coeffs_.c0_coeffs_ * car_state_.z1_memory_) +
alexbrandmeyer@668:    (car_coeffs_.a0_coeffs_ * car_state_.z2_memory_));
alexbrandmeyer@668:   car_state_.zy_memory_ = car_coeffs_.h_coeffs_ * car_state_.z2_memory_;
alexbrandmeyer@621:   // This section ripples the input-output path, to avoid delay...
alexbrandmeyer@621:   // It's the only part that doesn't get computed "in parallel":
alexbrandmeyer@668:   FPType in_out = input;
alexbrandmeyer@621:   for (int ch = 0; ch < n_ch_; ch++) {
alexbrandmeyer@620:     z1(ch) = z1(ch) + in_out;
alexbrandmeyer@621:     // This performs the ripple, and saves the final channel outputs in zy.
alexbrandmeyer@668:     in_out = car_state_.g_memory_(ch) * (in_out + car_state_.zy_memory_(ch));
alexbrandmeyer@668:     car_state_.zy_memory_(ch) = in_out;
alexbrandmeyer@620:   }
alexbrandmeyer@620:   car_state_.z1_memory_ = z1;
alexbrandmeyer@668:   *car_out = car_state_.zy_memory_;
alexbrandmeyer@620: }
alexbrandmeyer@620: 
alexbrandmeyer@621: // We start with a quadratic nonlinear function, and limit it via a
alexbrandmeyer@621: // rational function. This makes the result go to zero at high
alexbrandmeyer@620: // absolute velocities, so it will do nothing there.
alexbrandmeyer@668: void Ear::OHCNonlinearFunction(const FloatArray& velocities,
alexbrandmeyer@668:                                FloatArray* nonlinear_fun) {
alexbrandmeyer@668:   *nonlinear_fun = (1 + ((velocities * car_coeffs_.velocity_scale_) +
alexbrandmeyer@668:                          car_coeffs_.v_offset_).square()).inverse();
alexbrandmeyer@620: }
alexbrandmeyer@620: 
alexbrandmeyer@621: // This step is a one sample-time update of the inner-hair-cell (IHC) model,
alexbrandmeyer@621: // including the detection nonlinearity and either one or two capacitor state
alexbrandmeyer@621: // variables.
alexbrandmeyer@668: void Ear::IHCStep(const FloatArray& car_out, FloatArray* ihc_out) {
alexbrandmeyer@668:   FloatArray ac_diff = car_out - ihc_state_.ac_coupler_;
alexbrandmeyer@620:   ihc_state_.ac_coupler_ = ihc_state_.ac_coupler_ +
alexbrandmeyer@668:   (ihc_coeffs_.ac_coeff_ * ac_diff);
alexbrandmeyer@620:   if (ihc_coeffs_.just_hwr_) {
alexbrandmeyer@668:     FloatArray output(n_ch_);
alexbrandmeyer@668:     for (int ch = 0; ch < n_ch_; ++ch) {
alexbrandmeyer@668:       FPType a = (ac_diff(ch) > 0.0) ? ac_diff(ch) : 0.0;
alexbrandmeyer@668:       output(ch) = (a < 2) ? a : 2;
alexbrandmeyer@620:     }
alexbrandmeyer@668:     *ihc_out = output;
alexbrandmeyer@620:   } else {
alexbrandmeyer@668:     FloatArray conductance = CARFACDetect(ac_diff);
alexbrandmeyer@620:     if (ihc_coeffs_.one_cap_) {
alexbrandmeyer@668:       *ihc_out = conductance * ihc_state_.cap1_voltage_;
alexbrandmeyer@620:       ihc_state_.cap1_voltage_ = ihc_state_.cap1_voltage_ -
alexbrandmeyer@668:       (*ihc_out * ihc_coeffs_.out1_rate_) +
alexbrandmeyer@668:       ((1 - ihc_state_.cap1_voltage_)
alexbrandmeyer@668:        * ihc_coeffs_.in1_rate_);
alexbrandmeyer@620:     } else {
alexbrandmeyer@668:       *ihc_out = conductance * ihc_state_.cap2_voltage_;
alexbrandmeyer@620:       ihc_state_.cap1_voltage_ = ihc_state_.cap1_voltage_ -
alexbrandmeyer@668:       ((ihc_state_.cap1_voltage_ - ihc_state_.cap2_voltage_)
alexbrandmeyer@668:        * ihc_coeffs_.out1_rate_) +
alexbrandmeyer@668:       ((1 - ihc_state_.cap1_voltage_) * ihc_coeffs_.in1_rate_);
alexbrandmeyer@620:       ihc_state_.cap2_voltage_ = ihc_state_.cap2_voltage_ -
alexbrandmeyer@668:       (*ihc_out * ihc_coeffs_.out2_rate_) +
alexbrandmeyer@668:       ((ihc_state_.cap1_voltage_ - ihc_state_.cap2_voltage_)
alexbrandmeyer@668:        * ihc_coeffs_.in2_rate_);
alexbrandmeyer@620:     }
alexbrandmeyer@621:     // Here we smooth the output twice using a LPF.
alexbrandmeyer@668:     *ihc_out *= ihc_coeffs_.output_gain_;
alexbrandmeyer@620:     ihc_state_.lpf1_state_ = ihc_state_.lpf1_state_ +
alexbrandmeyer@668:     (ihc_coeffs_.lpf_coeff_ *
alexbrandmeyer@668:      (*ihc_out - ihc_state_.lpf1_state_));
alexbrandmeyer@620:     ihc_state_.lpf2_state_ = ihc_state_.lpf2_state_ +
alexbrandmeyer@668:     (ihc_coeffs_.lpf_coeff_ *
alexbrandmeyer@668:      (ihc_state_.lpf1_state_ - ihc_state_.lpf2_state_));
alexbrandmeyer@668:     *ihc_out = ihc_state_.lpf2_state_ - ihc_coeffs_.rest_output_;
alexbrandmeyer@620:   }
alexbrandmeyer@668:   ihc_state_.ihc_accum_ += *ihc_out;
alexbrandmeyer@620: }
alexbrandmeyer@620: 
alexbrandmeyer@668: bool Ear::AGCStep(const FloatArray& ihc_out) {
alexbrandmeyer@620:   int stage = 0;
alexbrandmeyer@668:   int n_stages = agc_coeffs_[0].n_agc_stages_;
alexbrandmeyer@668:   FPType detect_scale = agc_coeffs_[n_stages - 1].detect_scale_;
alexbrandmeyer@668:   bool updated = AGCRecurse(stage, detect_scale * ihc_out);
alexbrandmeyer@620:   return updated;
alexbrandmeyer@620: }
alexbrandmeyer@620: 
alexbrandmeyer@668: bool Ear::AGCRecurse(const int stage, FloatArray agc_in) {
alexbrandmeyer@620:   bool updated = true;
alexbrandmeyer@621:   // This is the decim factor for this stage, relative to input or prev. stage:
alexbrandmeyer@668:   int decim = agc_coeffs_[stage].decimation_;
alexbrandmeyer@621:   // This is the decim phase of this stage (do work on phase 0 only):
alexbrandmeyer@668:   int decim_phase = agc_state_[stage].decim_phase_ + 1;
alexbrandmeyer@620:   decim_phase = decim_phase % decim;
alexbrandmeyer@668:   agc_state_[stage].decim_phase_ = decim_phase;
alexbrandmeyer@621:   // Here we accumulate input for this stage from the previous stage:
alexbrandmeyer@668:   agc_state_[stage].input_accum_ += agc_in;
alexbrandmeyer@621:   // We don't do anything if it's not the right decim_phase.
alexbrandmeyer@621:   if (decim_phase == 0) {
alexbrandmeyer@621:     // Now we do lots of work, at the decimated rate.
alexbrandmeyer@621:     // These are the decimated inputs for this stage, which will be further
alexbrandmeyer@621:     // decimated at the next stage.
alexbrandmeyer@668:     agc_in = agc_state_[stage].input_accum_ / decim;
alexbrandmeyer@621:     // This resets the accumulator.
alexbrandmeyer@668:     agc_state_[stage].input_accum_ = FloatArray::Zero(n_ch_);
alexbrandmeyer@668:     if (stage < (agc_coeffs_.size() - 1)) {
alexbrandmeyer@621:       // Now we recurse to evaluate the next stage(s).
alexbrandmeyer@668:       // TODO (alexbrandmeyer): the Matlab version of AGCRecurse can return a
alexbrandmeyer@668:       // value for updated, but isn't used in that version. Check with Dick to
alexbrandmeyer@668:       // see if that is needed.
alexbrandmeyer@621:       updated = AGCRecurse(stage + 1, agc_in);
alexbrandmeyer@621:       // Afterwards we add its output to this stage input, whether it updated or
alexbrandmeyer@621:       // not.
alexbrandmeyer@668:       agc_in += agc_coeffs_[stage].agc_stage_gain_ *
alexbrandmeyer@668:       agc_state_[stage + 1].agc_memory_;
alexbrandmeyer@620:     }
alexbrandmeyer@668:     FloatArray agc_stage_state = agc_state_[stage].agc_memory_;
alexbrandmeyer@621:     // This performs a first-order recursive smoothing filter update, in time.
alexbrandmeyer@668:     agc_stage_state += agc_coeffs_[stage].agc_epsilon_ *
alexbrandmeyer@668:     (agc_in - agc_stage_state);
alexbrandmeyer@620:     agc_stage_state = AGCSpatialSmooth(stage, agc_stage_state);
alexbrandmeyer@668:     agc_state_[stage].agc_memory_ = agc_stage_state;
alexbrandmeyer@668:     updated = true;
alexbrandmeyer@620:   } else {
alexbrandmeyer@620:     updated = false;
alexbrandmeyer@620:   }
alexbrandmeyer@620:   return updated;
alexbrandmeyer@620: }
alexbrandmeyer@620: 
alexbrandmeyer@668: // TODO (alexbrandmeyer): figure out how to operate directly on stage_state.
alexbrandmeyer@668: // Using a pointer breaks the () indexing of the Eigen arrays, but there must
alexbrandmeyer@668: // be a way around this.
alexbrandmeyer@668: FloatArray Ear::AGCSpatialSmooth(const int stage, FloatArray stage_state) {
alexbrandmeyer@668:   int n_iterations = agc_coeffs_[stage].agc_spatial_iterations_;
alexbrandmeyer@620:   bool use_fir;
alexbrandmeyer@621:   use_fir = (n_iterations < 4) ? true : false;
alexbrandmeyer@620:   if (use_fir) {
alexbrandmeyer@668:     std::vector<FPType> fir_coeffs = agc_coeffs_[stage].agc_spatial_fir_;
alexbrandmeyer@620:     FloatArray ss_tap1(n_ch_);
alexbrandmeyer@620:     FloatArray ss_tap2(n_ch_);
alexbrandmeyer@620:     FloatArray ss_tap3(n_ch_);
alexbrandmeyer@620:     FloatArray ss_tap4(n_ch_);
alexbrandmeyer@668:     int n_taps = agc_coeffs_[stage].agc_spatial_n_taps_;
alexbrandmeyer@621:     // This initializes the first two taps of stage state, which are used for
alexbrandmeyer@621:     // both possible cases.
alexbrandmeyer@620:     ss_tap1(0) = stage_state(0);
alexbrandmeyer@621:     ss_tap1.block(1, 0, n_ch_ - 1, 1) = stage_state.block(0, 0, n_ch_ - 1, 1);
alexbrandmeyer@621:     ss_tap2(n_ch_ - 1) = stage_state(n_ch_ - 1);
alexbrandmeyer@621:     ss_tap2.block(0, 0, n_ch_ - 1, 1) = stage_state.block(1, 0, n_ch_ - 1, 1);
alexbrandmeyer@620:     switch (n_taps) {
alexbrandmeyer@620:       case 3:
alexbrandmeyer@668:         stage_state = (fir_coeffs[0] * ss_tap1) +
alexbrandmeyer@668:         (fir_coeffs[1] * stage_state) +
alexbrandmeyer@668:         (fir_coeffs[2] * ss_tap2);
alexbrandmeyer@620:         break;
alexbrandmeyer@620:       case 5:
alexbrandmeyer@621:         // Now we initialize last two taps of stage state, which are only used
alexbrandmeyer@621:         // for the 5-tap case.
alexbrandmeyer@620:         ss_tap3(0) = stage_state(0);
alexbrandmeyer@620:         ss_tap3(1) = stage_state(1);
alexbrandmeyer@621:         ss_tap3.block(2, 0, n_ch_ - 2, 1) =
alexbrandmeyer@668:         stage_state.block(0, 0, n_ch_ - 2, 1);
alexbrandmeyer@621:         ss_tap4(n_ch_ - 2) = stage_state(n_ch_ - 1);
alexbrandmeyer@621:         ss_tap4(n_ch_ - 1) = stage_state(n_ch_ - 2);
alexbrandmeyer@621:         ss_tap4.block(0, 0, n_ch_ - 2, 1) =
alexbrandmeyer@668:         stage_state.block(2, 0, n_ch_ - 2, 1);
alexbrandmeyer@668:         stage_state = (fir_coeffs[0] * (ss_tap3 + ss_tap1)) +
alexbrandmeyer@668:         (fir_coeffs[1] * stage_state) +
alexbrandmeyer@668:         (fir_coeffs[2] * (ss_tap2 + ss_tap4));
alexbrandmeyer@620:         break;
alexbrandmeyer@620:       default:
ronw@670:         assert(true && "Bad n_taps in AGCSpatialSmooth; should be 3 or 5.");
alexbrandmeyer@622:         break;
alexbrandmeyer@620:     }
alexbrandmeyer@620:   } else {
alexbrandmeyer@621:     stage_state = AGCSmoothDoubleExponential(stage_state,
alexbrandmeyer@668:                                              agc_coeffs_[stage].agc_pole_z1_,
alexbrandmeyer@668:                                              agc_coeffs_[stage].agc_pole_z2_);
alexbrandmeyer@620:   }
alexbrandmeyer@620:   return stage_state;
alexbrandmeyer@620: }
alexbrandmeyer@620: 
alexbrandmeyer@668: // TODO (alexbrandmeyer): figure out how to operate directly on stage_state.
alexbrandmeyer@668: // Same point as above for AGCSpatialSmooth.
alexbrandmeyer@621: FloatArray Ear::AGCSmoothDoubleExponential(FloatArray stage_state,
alexbrandmeyer@668:                                            const FPType pole_z1,
alexbrandmeyer@668:                                            const FPType pole_z2) {
alexbrandmeyer@622:   int32_t n_pts = stage_state.size();
alexbrandmeyer@621:   FPType input;
alexbrandmeyer@622:   FPType state = 0.0;
alexbrandmeyer@668:   // TODO (alexbrandmeyer): I'm assuming one dimensional input for now, but this
alexbrandmeyer@621:   // should be verified with Dick for the final version
alexbrandmeyer@668:   for (int i = n_pts - 11; i < n_pts; ++i){
alexbrandmeyer@621:     input = stage_state(i);
alexbrandmeyer@621:     state = state + (1 - pole_z1) * (input - state);
alexbrandmeyer@621:   }
alexbrandmeyer@668:   for (int i = n_pts - 1; i > -1; --i){
alexbrandmeyer@621:     input = stage_state(i);
alexbrandmeyer@621:     state = state + (1 - pole_z2) * (input - state);
alexbrandmeyer@621:   }
alexbrandmeyer@668:   for (int i = 0; i < n_pts; ++i){
alexbrandmeyer@621:     input = stage_state(i);
alexbrandmeyer@621:     state = state + (1 - pole_z1) * (input - state);
alexbrandmeyer@621:     stage_state(i) = state;
alexbrandmeyer@621:   }
alexbrandmeyer@620:   return stage_state;
alexbrandmeyer@620: }
alexbrandmeyer@621: 
alexbrandmeyer@668: FloatArray Ear::StageGValue(const FloatArray& undamping) {
alexbrandmeyer@668:   FloatArray r = car_coeffs_.r1_coeffs_ + car_coeffs_.zr_coeffs_ * undamping;
alexbrandmeyer@668:   return (1 - 2 * r * car_coeffs_.a0_coeffs_ + (r * r)) /
alexbrandmeyer@668:   (1 - 2 * r * car_coeffs_.a0_coeffs_ +
alexbrandmeyer@668:    car_coeffs_.h_coeffs_ * r * car_coeffs_.c0_coeffs_ + (r * r));
ronw@670: }