/home/marsyas/marsyas/src/marsyas/SeneffEar.cpp

/*
** Copyright (C) 1998-2006 George Tzanetakis <gtzan@cs.uvic.ca>
**  
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License as published by
** the Free Software Foundation; either version 2 of the License, or
** (at your option) any later version.
** 
** This program is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
** GNU General Public License for more details.
** 
** You should have received a copy of the GNU General Public License
** along with this program; if not, write to the Free Software 
** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/

#include "SeneffEar.h"
#include "fft.h"

#include <sstream>

using namespace std;
using namespace Marsyas;

SeneffEar::SeneffEar(mrs_string name):MarSystem("SeneffEar",name)
{
  //type_ = "SeneffEar";
  //name_ = name;

  firstUpdate = true;
  fs = 0.0f;
  stage = 4;

    addControls();
}

SeneffEar::~SeneffEar()
{
}

MarSystem* 
SeneffEar::clone() const 
{
  return new SeneffEar(*this);
}

void 
SeneffEar::addControls()
{
  addctrl("mrs_natural/stage", stage);
}

void 
SeneffEar::myUpdate(MarControlPtr sender)
{
    (void) sender;
  MRSDIAG("SeneffEar.cpp - SeneffEar:myUpdate");

  ostringstream name;

  //first time update
  if (firstUpdate) {
    //PreemphasisRTheta coefficients initialization
    stringstream matrix;
    matrix << "# MARSYAS realvec" << endl;
    matrix << "# Size = 8" << endl;
    matrix << endl;
    matrix << endl;
    matrix << "# type: matrix" << endl;
    matrix << "# rows: 4" << endl;
    matrix << "# columns: 2" << endl;
    matrix << "0.86 3.1148863";
    matrix << "0.99 0.0";
    matrix << "0.5  0.0";
    matrix << "0.95 3.14159";
    matrix << endl;
    matrix << "# Size = 8" << endl;
    matrix << "# MARSYAS realvec" << endl;
    matrix >> PreemphasisRThetaCoeffs;

    //FilterBankRTheta coefficients initialization
    matrix.clear();
    matrix << "# MARSYAS realvec" << endl;
    matrix << "# Size = 200" << endl;
    matrix << endl;
    matrix << endl;
    matrix << "# type: matrix" << endl;
    matrix << "# rows: 40" << endl;
    matrix << "# columns: 5" << endl;
    //         R-z       Theta-z       Radius       Theta       R-z2
    matrix << "0.0       3.14159       0.740055     2.633909    0.8";
    matrix << "0.86      2.997077      0.753637     2.178169    0.8";
    matrix << "0.86      2.879267      0.775569     1.856744    0.8";
    matrix << "0.86      2.761458      0.798336     1.617919    0.8";
    matrix << "0.86      2.643648      0.819169     1.433496    0.8";
    matrix << "0.86      2.525839      0.837158     1.286795    0.8";
    matrix << "0.8       2.964876      0.852598     1.167321    0.8";
    matrix << "0.86      2.408029      0.865429     1.068141    0.8";
    matrix << "0.86      2.29022       0.876208     0.984489    0.8";
    matrix << "0.86      2.17241       0.885329     0.912985    0.8";
    matrix << "0.86      2.054601      0.893116     0.851162    0.8";
    matrix << "0.86      1.936791      0.899823     0.797179    0.8";
    matrix << "0.8       2.788161      0.906118     0.749633    0.8";
    matrix << "0.86      1.818981      0.911236     0.70744     0.8";
    matrix << "0.86      1.701172      0.915747     0.669742    0.8";
    matrix << "0.86      1.583362      0.919753     0.635858    0.8";
    matrix << "0.86      1.465552      0.923335     0.605237    0.8";
    matrix << "0.86      1.347743      0.926565     0.57743     0.8";
    matrix << "0.8       2.611447      0.929914     0.552065    0.8";
    matrix << "0.86      1.229933      0.932576     0.528834    0.8";
    matrix << "0.86      1.112123      0.944589     0.487783    0.75";
    matrix << "0.86      0.994314      0.957206     0.452645    0.660714";
    matrix << "0.86      0.876504      0.956548     0.42223     0.672143";
    matrix << "0.86      0.758694      0.956653     0.395644    0.682143";
    matrix << "0.8       2.434732      0.956518     0.372208    0.690966";
    matrix << "0.86      0.640885      0.956676     0.351393    0.69881";
    matrix << "0.86      0.523075      0.956741     0.316044    0.712143";
    matrix << "0.8       2.258018      0.956481     0.287157    0.723052";
    matrix << "0.8       2.081304      0.956445     0.263108    0.732143";
    matrix << "0.8       1.904589      0.956481     0.242776    0.739835";
    matrix << "0.86      0.405265      0.958259     0.217558    0.749384";
    matrix << "0.8       1.727875      0.963083     0.197086    0.757143";
    matrix << "0.8       1.55116       0.969757     0.175115    0.769048";
    matrix << "0.8       1.374446      0.97003      0.153697    0.780662";
    matrix << "0.8       1.197732      0.970382     0.134026    0.791337";
    matrix << "0.8       1.021017      0.970721     0.118819    0.799596";
    matrix << "0.8       1.5           0.970985     0.106711    0.8";
    matrix << "0.8       1.2           0.971222     0.096843    0.8";
    matrix << "0.8       1.0           0.97144      0.088645    0.8";
    matrix << "0.8       0.9           0.971645     0.081727    0.8";
    matrix << endl;
    matrix << "# Size = 200" << endl;
    matrix << "# MARSYAS realvec" << endl;
    matrix >> FilterBankRThetaCoeffs;

    realvec a,b,c;

    //SeneffPreemphasis coefficients computation
    a.create(3);
    b.create(3);
    a(0) = PreemphasisRThetaCoeffs(0,0)*PreemphasisRThetaCoeffs(0,0);
    a(1) = -2.0f*PreemphasisRThetaCoeffs(0,0)*cos(PreemphasisRThetaCoeffs(0,1));
    a(2) = 1.0f;
    b(0) = PreemphasisRThetaCoeffs(1,0)*PreemphasisRThetaCoeffs(1,0);
    b(1) = -2.0f*PreemphasisRThetaCoeffs(1,0)*cos(PreemphasisRThetaCoeffs(1,1));
    b(2) = 1.0f;
    polyConv(a, b, a);
    b(0) = PreemphasisRThetaCoeffs(2,0)*PreemphasisRThetaCoeffs(2,0);
    b(1) = -2.0f*PreemphasisRThetaCoeffs(2,0)*cos(PreemphasisRThetaCoeffs(2,1));
    b(2) = 1.0f;
    polyConv(a, b, a);
    b(0) = PreemphasisRThetaCoeffs(3,0)*PreemphasisRThetaCoeffs(3,0);
    b(1) = -2.0f*PreemphasisRThetaCoeffs(3,0)*cos(PreemphasisRThetaCoeffs(3,1));
    b(2) = 1.0f;
    polyConv(a, b, c);
    SeneffPreemphasisCoeffs.create(c.getSize());
      
    for (mrs_natural i = 0; i < c.getSize(); ++i)
        SeneffPreemphasisCoeffs(i) = c(c.getSize() - i - 1);

    //SeneffFilterBank coefficients computation
    channels = FilterBankRThetaCoeffs.getRows();
    SeneffFilterBankCoeffs.create(channels, 5);
    mrs_real sum;
    for (mrs_natural i = 0; i < channels; ++i) {
      SeneffFilterBankCoeffs(i,0) = 1.0f;
      SeneffFilterBankCoeffs(i,1) = -2.0f*FilterBankRThetaCoeffs(i,0)*cos(FilterBankRThetaCoeffs(i,1));
      SeneffFilterBankCoeffs(i,2) = FilterBankRThetaCoeffs(i,0)*FilterBankRThetaCoeffs(i,0);
      sum = SeneffFilterBankCoeffs(i,0) + SeneffFilterBankCoeffs(i,1) + SeneffFilterBankCoeffs(i,2);
      SeneffFilterBankCoeffs(i,0) /= sum;
      SeneffFilterBankCoeffs(i,1) /= sum;
      SeneffFilterBankCoeffs(i,2) /= sum;
    }

    //SeneffForward & SeneffBackward coefficients computation
    SeneffForwardCoeffs.create(mrs_natural(5), static_cast<mrs_natural>(channels));
    SeneffBackwardCoeffs.create(mrs_natural(5), static_cast<mrs_natural>(channels));
    a.create(3),b.create(3);
    for (mrs_natural j = 0; j < channels; j++){
      a(0) = FilterBankRThetaCoeffs(j,4)*FilterBankRThetaCoeffs(j,4);
      a(1) = -2.0f*FilterBankRThetaCoeffs(j,4)*cos(FilterBankRThetaCoeffs(j,3)/2);
      a(2) = 1.0f;
      polyConv(a, a, c);
      for (mrs_natural i = 0; i < 5; ++i) SeneffForwardCoeffs(i,j) = c(5 - i - 1);
      b(0) = FilterBankRThetaCoeffs(j,2)*FilterBankRThetaCoeffs(j,2);
      b(1) = -2.0f*FilterBankRThetaCoeffs(j,2)*cos(FilterBankRThetaCoeffs(j,3));
      b(2) = 1.0f;
      polyConv(b, b, c);
      for (mrs_natural i = 0; i < 5; ++i) SeneffBackwardCoeffs(i,j) = c(5 - i - 1);
    }

    //SeneffPreemphasis filter creation and coefficients initialization
    SeneffPreemphasisFilter = new Filter("SeneffPreemphasisFilter");
    a.create(1);
    a(0) = 1.0f;
    SeneffPreemphasisFilter->setctrl("mrs_realvec/ncoeffs", SeneffPreemphasisCoeffs);
    SeneffPreemphasisFilter->setctrl("mrs_realvec/dcoeffs", a);

    //SeneffFilterBank filters cascade creation and coefficients initialization
    SeneffFilterBank = new Cascade("SeneffFilterBank");
    Filter* filter;
    a.create(3);
    b.create(3);
    for (mrs_natural i = 0; i < channels; ++i) {
      name.clear();
      name.str("");
      name << "filter_" << i;
      filter = new Filter(name.str());
      b(0) = SeneffFilterBankCoeffs(i,0);
      b(1) = SeneffFilterBankCoeffs(i,1);
      b(2) = SeneffFilterBankCoeffs(i,2);
      a(0) = 1.0f;
      a(1) = SeneffFilterBankCoeffs(i,3);
      a(2) = SeneffFilterBankCoeffs(i,4);
      filter->setctrl("mrs_realvec/ncoeffs", b);
      filter->setctrl("mrs_realvec/dcoeffs", a);
      SeneffFilterBank->addMarSystem(filter);
    }

    //resonatorFilter creation
    resonatorFilter = new Parallel("resonatorFilter");
    a.create(5);
    b.create(5);
    for (mrs_natural i = 0; i < channels; ++i) {
      name.clear();
      name.str("");
      name << "filter_" << i;
      filter = new Filter(name.str());
      for (mrs_natural j = 0; j < 5; j++) b(j) = SeneffForwardCoeffs(j,i);
      for (mrs_natural j = 0; j < 5; j++) a(j) = SeneffBackwardCoeffs(j,i);
      filter->setctrl("mrs_realvec/ncoeffs", b);
      filter->setctrl("mrs_realvec/dcoeffs", a);
      resonatorFilter->addMarSystem(filter);
    }

    //run an impulse through the preemphasis filter
    realvec impulse;
    impulse.create(mrs_natural(1), mrs_natural(256));
    impulse(0) = 1.0f;
    realvec y0;
    y0.create(impulse.getRows(), impulse.getCols());
    SeneffPreemphasisFilter->setctrl("mrs_natural/inSamples", impulse.getCols());
    SeneffPreemphasisFilter->setctrl("mrs_natural/inObservations", impulse.getRows());
    SeneffPreemphasisFilter->update();
    SeneffPreemphasisFilter->process(impulse, y0);

    //then through the cascade of zeros
    realvec y1;
    y1.create(y0.getRows()*channels, y0.getCols());
    SeneffFilterBank->setctrl("Filter/filter_0/mrs_natural/inSamples", y0.getCols());
    SeneffFilterBank->setctrl("Filter/filter_0/mrs_natural/inObservations", y0.getRows());
    SeneffFilterBank->update();
    SeneffFilterBank->process(y0, y1);

    //Finally through each of the resonator filters
    y.create(y1.getRows(), y1.getCols());
    resonatorFilter->setctrl("Filter/filter_0/mrs_natural/inSamples", y1.getCols());
    for (mrs_natural i = 0; i < channels; ++i) {
      name.clear();
      name.str("");
      name << "Filter/" << "filter_" << i << "/mrs_natural/inObservations";
      resonatorFilter->setctrl(name.str(), y0.getRows());
    }
    resonatorFilter->update();
    resonatorFilter->process(y1, y);

    //reset the state vector of each filter
    realvec state;
    state.create(SeneffPreemphasisFilter->getctrl("mrs_realvec/state")->to<mrs_realvec>().getRows(),SeneffPreemphasisFilter->getctrl("mrs_realvec/state")->to<mrs_realvec>().getCols());
    SeneffPreemphasisFilter->setctrl("mrs_natural/stateUpdate", mrs_natural(1));
    SeneffPreemphasisFilter->updControl("mrs_realvec/state", state);
    SeneffPreemphasisFilter->setctrl("mrs_natural/stateUpdate", mrs_natural(0));

    state.create(SeneffFilterBank->getctrl("Filter/filter_0/mrs_realvec/state")->to<mrs_realvec>().getRows(),SeneffFilterBank->getctrl("Filter/filter_0/mrs_realvec/state")->to<mrs_realvec>().getCols());
    for (mrs_natural i = 0; i < channels; ++i) {
      name.clear();
      name.str("");
      name << "Filter/" << "filter_" << i << "/mrs_natural/stateUpdate";
      SeneffFilterBank->setctrl(name.str(), mrs_natural(1));
      name.clear();
      name.str("");
      name << "Filter/" << "filter_" << i << "/mrs_realvec/state";
      SeneffFilterBank->updControl(name.str(), state);
      name.clear();
      name.str("");
      name << "Filter/" << "filter_" << i << "/mrs_natural/stateUpdate";
      SeneffFilterBank->setctrl(name.str(), mrs_natural(0));
    }

    state.create(resonatorFilter->getctrl("Filter/filter_0/mrs_realvec/state")->to<mrs_realvec>().getRows(),resonatorFilter->getctrl("Filter/filter_0/mrs_realvec/state")->to<mrs_realvec>().getCols());
    for (mrs_natural i = 0; i < channels; ++i) {
      name.clear();
      name.str("");
      name << "Filter/" << "filter_" << i << "/mrs_natural/stateUpdate";
      resonatorFilter->setctrl(name.str(), mrs_natural(1));
      name.clear();
      name.str("");
      name << "Filter/" << "filter_" << i << "/mrs_realvec/state";
      resonatorFilter->updControl(name.str(), state);
      name.clear();
      name.str("");
      name << "Filter/" << "filter_" << i << "/mrs_natural/stateUpdate";
      resonatorFilter->setctrl(name.str(), mrs_natural(0));
    }

    //constants initialization
    hwrA = 10.0f;
    hwrB = 65.0f;
    hwrG = 2.35f;

    Cn.create(channels);

    //lowPassFilter creation and coefficients initialization
    lpAlpha = 0.209611f;
    lowPassFilter = new Filter("lowPassFilter");
    a.create(2);
    a(0) = -lpAlpha;
    a(1) = 1.0f;
    polyConv(a, a, a);
    polyConv(a, a, a);
    polyFlip(a);
    b.create(1);
    b(0) = a.sum();
    lowPassFilter->setctrl("mrs_realvec/ncoeffs", b);
    lowPassFilter->setctrl("mrs_realvec/dcoeffs", a);

    //Seneff's Adaptive Gain Control filter creation and coefficients initialization
    initial_yn = (mrs_real)0.23071276;
    alpha_agc = (mrs_real)0.979382181;
    kagc = (mrs_real)0.002;
    AGCfilter = new Filter("AGCfilter");
    a.create(1);
    a(0) = alpha_agc;
    b.create(2);
    b(1) = 1.0f - alpha_agc;
    AGCfilter->setctrl("mrs_realvec/ncoeffs", b);
    AGCfilter->setctrl("mrs_realvec/dcoeffs", a);
    state.create(channels, mrs_natural(1));
    state.setval(initial_yn);
    AGCfilter->setctrl("mrs_natural/inObservations", channels);
    AGCfilter->setctrl("mrs_natural/stateUpdate", mrs_natural(1));
    AGCfilter->updControl("mrs_realvec/state", state);
    AGCfilter->setctrl("mrs_natural/stateUpdate", mrs_natural(0));

    firstUpdate = false;
  }//end first update

  //--------------------------------------------------------------------------------
  //compute the FFT then find the gain peak then divide each forward polynomial by the maximum gain (to normalize) and then multiply by the desired low frequency roll-off.
  if (fs != getctrl("mrs_real/israte")->to<mrs_real>()) {
    fs = getctrl("mrs_real/israte")->to<mrs_real>();
    fft fft;
    realvec Y(y.getCols());
    mrs_real *Yd = Y.getData();
    mrs_real abs2(1.0f),maxAbs2(0.0f);
    mrs_real gain;
    mrs_real rolloff;
    SeneffForwardCoeffsNormalized.create(SeneffForwardCoeffs.getRows(), SeneffForwardCoeffs.getCols());
    for (mrs_natural i = 0; i < channels; ++i) {
      for (mrs_natural j = 0; j < Y.getCols(); j++) Y(j) = y(i,j);
      fft.rfft(Yd, Y.getCols()/2, FFT_FORWARD);
      maxAbs2 = Y(0)*Y(0);
      for (mrs_natural j = 1; j < y.getCols()/2; j++) if ((abs2 = Y(2*j)*Y(2*j) + Y(2*j+1)*Y(2*j+1)) > maxAbs2) maxAbs2 = abs2;
      gain = 1/(y.getCols()*sqrt(maxAbs2));
      rolloff = min((mrs_real)((FilterBankRThetaCoeffs(i,3)/PI*fs/2)/1600),(mrs_real)1.0);
      for (mrs_natural j = 0; j < SeneffForwardCoeffsNormalized.getRows(); j++) SeneffForwardCoeffsNormalized(j,i) = SeneffForwardCoeffs(j,i)*gain*rolloff;
    }
    //then update the resonatorFilter
    realvec b(5);
    for (mrs_natural i = 0; i < channels; ++i) {
      name.clear();
      name.str("");
      name << "Filter/" << "filter_" << i << "/mrs_realvec/ncoeffs";
      for (mrs_natural j = 0; j < 5; j++) b(j) = SeneffForwardCoeffsNormalized(j,i);
      resonatorFilter->setctrl(name.str(), b);
    }
    resonatorFilter->update();

    Tua = (mrs_real)(58.3333/fs);
    Tub = (mrs_real)(8.3333/fs);
  }

  //--------------------------------------------------------------------------------
  //SeneffPreemphasis update
  SeneffPreemphasisFilter->setctrl("mrs_natural/inSamples", getctrl("mrs_natural/inSamples"));
  SeneffPreemphasisFilter->setctrl("mrs_natural/inObservations", getctrl("mrs_natural/inObservations"));
  SeneffPreemphasisFilter->setctrl("mrs_real/israte", getctrl("mrs_real/israte"));
  SeneffPreemphasisFilter->update();
  if ((int)slice_0.getSize() != SeneffPreemphasisFilter->getctrl("mrs_natural/onObservations")->to<mrs_natural>() * SeneffPreemphasisFilter->getctrl("mrs_natural/onSamples")->to<mrs_natural>()) {
    slice_0.create(SeneffPreemphasisFilter->getctrl("mrs_natural/onObservations")->to<mrs_natural>(), SeneffPreemphasisFilter->getctrl("mrs_natural/onSamples")->to<mrs_natural>());
  }

  //SeneffFilterBank filters cascade update
  SeneffFilterBank->setctrl("Filter/filter_0/mrs_natural/inSamples", SeneffPreemphasisFilter->getctrl("mrs_natural/onSamples"));
  SeneffFilterBank->setctrl("Filter/filter_0/mrs_natural/inObservations", SeneffPreemphasisFilter->getctrl("mrs_natural/onObservations"));
  SeneffFilterBank->setctrl("Filter/filter_0/mrs_real/israte", SeneffPreemphasisFilter->getctrl("mrs_real/osrate"));
  SeneffFilterBank->update();
  if ((int)slice_1.getSize() != SeneffFilterBank->getctrl("mrs_natural/onObservations")->to<mrs_natural>() * SeneffFilterBank->getctrl("mrs_natural/onSamples")->to<mrs_natural>()) {
    slice_1.create(SeneffFilterBank->getctrl("mrs_natural/onObservations")->to<mrs_natural>(), SeneffFilterBank->getctrl("mrs_natural/onSamples")->to<mrs_natural>());
  }

  //resonatorFilter update
  resonatorFilter->setctrl("Filter/filter_0/mrs_natural/inSamples", SeneffFilterBank->getctrl("mrs_natural/onSamples"));
  for (mrs_natural i = 0; i < channels; ++i) {
    name.clear();
    name.str("");
    name << "Filter/" << "filter_" << i << "/mrs_natural/inObservations";
    resonatorFilter->setctrl(name.str(), SeneffFilterBank->getctrl("Filter/filter_0/mrs_natural/onObservations"));
  }
  resonatorFilter->setctrl("Filter/filter_0/mrs_real/israte", SeneffFilterBank->getctrl("mrs_real/osrate"));
  resonatorFilter->update();
  if ((int)slice_2.getSize() != resonatorFilter->getctrl("mrs_natural/onObservations")->to<mrs_natural>() * resonatorFilter->getctrl("mrs_natural/onSamples")->to<mrs_natural>()) {
    slice_2.create(resonatorFilter->getctrl("mrs_natural/onObservations")->to<mrs_natural>(), resonatorFilter->getctrl("mrs_natural/onSamples")->to<mrs_natural>());
  }

  //lowPassFilter update
  lowPassFilter->setctrl("mrs_natural/inSamples", resonatorFilter->getctrl("mrs_natural/onSamples"));
  lowPassFilter->setctrl("mrs_natural/inObservations", resonatorFilter->getctrl("mrs_natural/onObservations"));
  lowPassFilter->setctrl("mrs_real/israte", resonatorFilter->getctrl("mrs_real/osrate"));
  lowPassFilter->update();
  if ((int)slice_3.getSize() != lowPassFilter->getctrl("mrs_natural/onObservations")->to<mrs_natural>() * lowPassFilter->getctrl("mrs_natural/onSamples")->to<mrs_natural>()) {
    slice_3.create(lowPassFilter->getctrl("mrs_natural/onObservations")->to<mrs_natural>(), lowPassFilter->getctrl("mrs_natural/onSamples")->to<mrs_natural>());
  }

  //AGCfilter update
  AGCfilter->setctrl("mrs_natural/inSamples", lowPassFilter->getctrl("mrs_natural/onSamples"));
  AGCfilter->setctrl("mrs_natural/inObservations", lowPassFilter->getctrl("mrs_natural/onObservations"));
  AGCfilter->setctrl("mrs_real/israte", lowPassFilter->getctrl("mrs_real/osrate"));
  AGCfilter->update();

  //this update
  setctrl("mrs_natural/onSamples", AGCfilter->getctrl("mrs_natural/onSamples"));
  setctrl("mrs_natural/onObservations", AGCfilter->getctrl("mrs_natural/onObservations")->to<mrs_natural>());
  setctrl("mrs_real/osrate", AGCfilter->getctrl("mrs_real/osrate"));
}

void 
SeneffEar::polyConv(realvec& a, realvec& b, realvec& c)
{
  mrs_natural la(a.getSize());
  mrs_natural lb(b.getSize());
  mrs_natural n = la + lb - 1;

  realvec ta(a); ta.stretch(n);
  realvec tb(b); tb.stretch(n);
  realvec tc;    tc.create(n);

  for (mrs_natural k = 0; k < n; k++){
    for (mrs_natural i = 0; i <= k; ++i){
      tc(k) += ta(i)*tb(k-i);
    }
  }
  if ((mrs_natural)c.getSize() != n) c.create(n);
  c = tc;
}

void 
SeneffEar::polyFlip(realvec& a)
{
  mrs_natural la(a.getSize());
  realvec ta(a);

  for (mrs_natural i = 0; i < la; ++i){
    a(i) = ta(la - i - 1);
  }
}


void 
SeneffEar::myProcess(realvec& in, realvec& out)
{
  checkFlow(in, out);
  //lmartins: if (mute_) return;
    if(getctrl("mrs_bool/mute")->to<mrs_bool>()) return;

  mrs_natural s = 0;
  stage = getctrl("mrs_natural/stage")->to<mrs_natural>();

  SeneffPreemphasisFilter->process(in, slice_0);
  SeneffFilterBank->process(slice_0, slice_1);
  if (s++ == stage) {out = slice_1; return;}
  resonatorFilter->process(slice_1, slice_2);
  if (s++ == stage) {out = slice_2; return;}
  //Seneff's detector non-linearity
  for (mrs_natural i = 0; i < slice_2.getRows(); ++i) {
    for (mrs_natural j = 0; j < slice_2.getCols(); j++) {
      slice_2(i,j) = hwrA*atan(hwrB*max((mrs_real)0.0,slice_2(i,j))) + exp(hwrA*hwrB*min((mrs_real)0.0,slice_2(i,j)));
    }
  }
  //Seneff's short-term adaptation (a reservoir hair cell model)
  mrs_real flow;
  for (mrs_natural j = 0; j < slice_2.getCols(); j++) {
    for (mrs_natural i = 0; i < slice_2.getRows(); ++i) {
      flow = max((mrs_real)0.0,Tua*(slice_2(i,j)-Cn(i)));
      Cn(i) = Cn(i) + flow - Tub*Cn(i);
      slice_2(i,j) = flow;
    }
  }
  if (s++ == stage) {out = slice_2; return;}
  //Seneff's Low Pass filter (to model the loss of Synchrony)
  lowPassFilter->process(slice_2, slice_3);
  if (s++ == stage) {out = slice_3; return;}
  //Seneff's Adaptive Gain Control
  AGCfilter->process(slice_3, out);
  for (mrs_natural i = 0; i < out.getRows(); ++i) {
    for (mrs_natural j = 0; j < out.getCols(); j++) {
      out(i,j) = slice_3(i,j)/(1.0f + kagc*out(i,j));
    }
  }
}

---