// To do : 1. change the option that about # of model
// 2. add the option that can generate the synapse output / spike rate
// 3. check the parameter for CAT/HUMAN, and maybe more
// 4. Species can be specified as others as the same in the paper
/* Model numbers as used in ARLO Heinz et al., Fig. 4 */
// model == 1 bmmodel = FeedForward_NL;
// model == 2 bmmodel = FeedBack_NL;
// model == 3 bmmodel = Sharp_Linear;
// model == 4 bmmodel = Broad_Linear;
// model == 5 bmmodel = Broad_Linear_High;
#include "cmpa.h"
#include "hc.h"
#include "synapse.h"
#include "complex.h"
#include "filters.h"
#include <stdlib.h>
double Get_tau(int species, double cf, int order, double* taumax, double* taumin, double* taurange);
double erbGM(double);
double cochlea_f2x(int species,double f);
double cochlea_x2f(int species,double x);
void runAN2(TAuditoryNerve *p,const double *in,double *out,const int length);
double runBasilarMembrane(TBasilarMembrane *bm, double x);
void run2BasilarMembrane(TBasilarMembrane *bm, const double *in, double *out, const int length);
/*
* different species have different parameters !!!!
* 1. Get_Tau(...)
* 2. Synapse setup in this function
* 3. IHC low pass filter cutoff freq
*
* */
void initAuditoryNerve(TAuditoryNerve *p,int model, int species, double tdres, double cf, double spont)
{ /* default */
Tsynapse *syn;
double Kcf,temp;
double p1,p2,pst,psl;
p->run = runAN;
p->run2 = runAN2;
syn = &(p->syn);
p->model = model;
p->species = species;
p->cf = cf;
p->spont = spont;
p->tdres = tdres;
p->run = runAN;
p->run2 = runAN2;
/*
* Init the basilar membrane
* */
initBasilarMembrane(&(p->bm),model,species,tdres,cf);
/* Now add support for different species */
/* This line outputs the appropriate sout rate based on whether or not spikes are being generated (with refractoriness) */
if(p->ifspike==1) syn->Ass = 350;
else syn->Ass = 130; /* borrow this for Mike G. Heinz */
switch(species)
{
case 1 : /* Cat - low CFs only (Carney '93) */
case 9 : /* Cat - all CFs "Universal" (Zhang et al 2001) */
/*
* Init the synapse
* */
syn->tdres = tdres;
syn->cf = cf;
syn->spont = spont;
syn->Pimax = 0.6;
syn->PTS = 1+9*spont/(9+spont);
syn->Ar_over_Ast = 6.0;
syn->tauR = 0.002;
syn->tauST = 0.060;
initSynapse(syn); /*This function calcuate all other parameters as described in the appendix*/
/*
* Init the ihcppi rectify function and parameters
* !!!!Use the result from initSynapse(syn)
*
* */
Kcf = 2.0+3.0*log10(cf/1000);
if(Kcf<1.5) Kcf = 1.5;
syn->Vsat = syn->Pimax*Kcf*20.0*(1+spont)/(5+spont);
/* Important !!! The original code seems use the following equation!!!
This is not shown in the paper but Prest is very small compare to the other part < 0.01*Vsat
and the results is the same
*/
/* syn->Vsat = syn->Pimax*Kcf*20.0*(1+spont)/(5+spont)+syn->Prest; */
break;
case 0 :
default: /* Human */
/*
* Init the synapse
* */
syn->tdres = tdres;
syn->cf = cf;
syn->spont = spont;
syn->Pimax = 0.6;
syn->PTS = 1+9*spont/(9+spont); /* Peak to Steady State Ratio, characteristic of PSTH */
syn->Ar_over_Ast = 6.0;
syn->tauR = 0.002;
syn->tauST = 0.060;
initSynapse(syn); /*This function calcuate all other parameters as described in the appendix*/
/*
* Init the ihcppi rectify function and parameters
* !!!!Use the result from initSynapse(syn)
*
* */
Kcf = 2.0+1.3*log10(cf/1000); /*/ From M. G. Heinz */
if(Kcf<1.5) Kcf = 1.5;
syn->Vsat = syn->Pimax*Kcf*60.0*(1+spont)/(6+spont);
/* Important !!! The original code seems use the following equation!!!
This is not shown in the paper but Prest is very small compare to the other part < 0.01*Vsat
and the results is the same, so we use the upper equation
*/
/* syn->Vsat = syn->Pimax*Kcf*60.0*(1+spont)/(6+spont)+syn->Prest; */
break;
};
temp = log(2)*syn->Vsat/syn->Prest;
if(temp<400) pst = log(exp(temp)-1);
else pst = temp;
psl = syn->Prest*pst/log(2);
p2 = pst;
p1 = psl/pst;
#ifdef DEBUG
printf("\ninitSyanpse : p1=%f; p2=%f",p1,p2);
#endif
p->ihcppi.psl = psl;
p->ihcppi.pst = pst;
p->ihcppi.p1 = p1;
p->ihcppi.p2 = p2;
p->ihcppi.run = runIHCPPI; /*These are functions used to get ihcppi output */
p->ihcppi.run2 = runIHCPPI2; /*Defined in hc.c, hc.h */
/*
* Init Inner Hair Cell Model
*/
/* For Human Only !!!!!!! Hair Cell low pass filter (tdres,cutoff,gain,order) */
switch(species)
{
case 0:
default: /* Human */
initLowPass(&(p->ihc.hclp),tdres,4500.0,1.0,7);
break;
case 1: /* Cat, this version use the Laurel's revcor data to determine taumin, taumax - only good for low freqs*/
case 9: /* Also for universal CAT; as in JASA 2001 Zhang et al. - good for all CFs */
initLowPass(&(p->ihc.hclp),tdres,3800.0,1.0,7);
break;
};
p->ihc.run = runHairCell;
p->ihc.run2 = runHairCell2;
p->ihc.hcnl.A0 = 0.1; /* Inner Hair Cell Nonlinear Function */
p->ihc.hcnl.B = 2000.;
p->ihc.hcnl.C = 1.74;
p->ihc.hcnl.D = 6.87e-9;
p->ihc.hcnl.run = runIHCNL;
p->ihc.hcnl.run2 = runIHCNL2;
return;
}
double runAN(TAuditoryNerve *p,double x)
{
double bmout,ihcout,ppiout,sout;
bmout = p->bm.run(&(p->bm),x);
ihcout = p->ihc.run(&(p->ihc),bmout);
ppiout = p->ihcppi.run(&(p->ihcppi),ihcout);
sout = p->syn.run(&(p->syn),ppiout);
return(sout);
}
void runAN2(TAuditoryNerve *p,const double *in,double *out,const int length)
{
p->bm.run2(&(p->bm),in,out,length);
p->ihc.run2(&(p->ihc),out,out,length);
p->ihcppi.run2(&(p->ihcppi),out,out,length);
p->syn.run2(&(p->syn),out,out,length);
return;
}
/* User Get_Tau, initGammaTone, initBoltzman*/
void initBasilarMembrane(TBasilarMembrane* bm,int model, int species, double tdres, double cf)
{ /*
*
* ##### Get Basilar Membrane ########
* 1. Get a structure of BasilarMembrane
* 2. Specify wide band filter in the BasilarMembrane Model
* //WB filter not used for all model versions, but it's computed here anyway
* 3. Specify the OHC model in the control path
* 3.2 Specify the NL function used in the OHC
* 4. Specify the AfterOHC NL in the control path
* */
int bmmodel;
double taumax,taumin,taurange; /* general */
double x, centerfreq,tauwb,tauwbmin,dtmp,wb_gain; /* for wb */
double absdb,ohcasym; /* for ohc */
double dc,R,minR; /* for afterohc */
TNonLinear *p;
/* Model numbers as used in ARLO Heinz et al., Fig. 4 */
if(model == 1) bmmodel = FeedForward_NL;
else if(model == 2) bmmodel = FeedBack_NL;
else if(model == 3) bmmodel = Sharp_Linear;
else if(model == 4) bmmodel = Broad_Linear;
else if(model == 5) bmmodel = Broad_Linear_High;
bm->bmmodel = bmmodel;
bm->tdres = tdres;
/*
* Determine taumax,taumin,order here
*/
bm->run = runBasilarMembrane;
bm->run2 = run2BasilarMembrane;
bm->bmorder = 3;
Get_tau(species,cf,bm->bmorder,&taumax,&taumin,&taurange);
bm->TauMax = taumax;
bm->TauMin = taumin;
if(bm->bmmodel&Broad_ALL)
bm->tau = taumin;
else
bm->tau = taumax;
initGammaTone(&(bm->bmfilter),tdres,cf,bm->tau,1.0,bm->bmorder);
initGammaTone(&(bm->gfagain),tdres,cf,taumin,1.0,1);
/*
* Get Wbfilter parameters
*/
x = cochlea_f2x(species,cf);
centerfreq = cochlea_x2f(species,x+1.2); /* shift the center freq Qing use 1.1 shift */
bm->wborder = 3;
tauwb = taumin+0.2*(taumax-taumin);
tauwbmin = tauwb/taumax*taumin;
dtmp = tauwb*TWOPI*(centerfreq-cf);
wb_gain = pow((1+dtmp*dtmp), bm->wborder/2.0);
bm->TauWB = tauwb;
bm->TauWBMin = tauwbmin;
initGammaTone(&(bm->wbfilter), tdres, centerfreq, tauwb, wb_gain, bm->wborder);
bm->A=(tauwb/taumax-tauwbmin/taumin)/(taumax-taumin);
bm->B=(taumin*taumin*tauwb-taumax*taumax*tauwbmin)/(taumax*taumin*(taumax-taumin));
/*
* Init OHC model
*/
absdb = 20; /* The value that the BM starts compression */
initLowPass(&(bm->ohc.hclp), tdres, 800.0, 1.0, 3); /* This is now use in both Human & Cat MODEL */
/*/ parameter into boltzman is corner,slope,strength,x0,s0,x1,s1,asym
// The corner determines the level that BM have compressive nonlinearity */
ohcasym = 7.0;
/*/set OHC Nonlinearity as boltzman function combination */
init_boltzman(&(bm->ohc.hcnl),absdb-12,0.22,0.08,5,12,5,5,ohcasym);
bm->ohc.run = runHairCell;
bm->ohc.run2 = runHairCell2;
/*
* Init AfterOHC model
*/
p = &(bm->afterohc);
dc = (ohcasym-1)/(ohcasym+1.0)/2.0;
dc-=0.05;
minR = 0.05;
p->TauMax = taumax;
p->TauMin = taumin;
R = taumin/taumax;
if(R<minR) p->minR = 0.5*R;
else p->minR = minR;
p->A = p->minR/(1-p->minR); /* makes x = 0; output = 1; */
p->dc = dc;
R = R-p->minR;
/* This is for new nonlinearity */
p->s0 = -dc/log(R/(1-p->minR));
p->run = runAfterOhcNL;
p->run2 = runAfterOhcNL2;
return;
}
/** pass the signal through the tuning filter.
using different model
Sharp_Linear | Broad_Linear | Broad_Linear_High | FeedBack_NL | FeedForward_NL
*/
double runBasilarMembrane(TBasilarMembrane *bm, double x){
double out;
const int length = 1;
run2BasilarMembrane(bm, &x, &out, length);
return(out);
}
void run2BasilarMembrane(TBasilarMembrane *bm, const double *in, double *out, const int length)
{
register int i;
double wb_gain,dtmp,taunow;
double wbout,ohcout;
double x, x1,out1;
for(i=0; i<length; i++)
{
x = in[i];
x1 = bm->bmfilter.run(&(bm->bmfilter),x); /*/ pass the signal through the tuning filter */
switch(bm->bmmodel){
default:
case Sharp_Linear: /* / if linear model, needn't get the control signal */
case Broad_Linear:
out1 = x1; break;
case Broad_Linear_High: /* /adjust the gain of the tuning filter */
out1 = x1*pow((bm->TauMin/bm->TauMax),bm->bmfilter.Order);
break;
case FeedBack_NL: /* /get the output of the tuning filter as the control signal */
wbout = x1;
ohcout = bm->ohc.run(&(bm->ohc),wbout); /*/ pass the control signal through OHC model */
/*/ pass the control signal through nonliearity after OHC */
bm->tau = bm->afterohc.run(&(bm->afterohc),ohcout);
/*/ set the tau of the tuning filter */
bm->bmfilter.settau(&(bm->bmfilter),bm->tau);
/* Gain Control of the tuning filter */
out1 = pow((bm->tau/bm->TauMax),bm->bmfilter.Order)*x1;
break;
case FeedForward_NL:
/*/get the output of the wide-band pass as the control signal */
wbout = bm->wbfilter.run(&(bm->wbfilter),x);
/*/ scale the tau for wide band filter in control path */
taunow = bm->A*bm->tau*bm->tau-bm->B*bm->tau;
bm->wbfilter.settau(&(bm->wbfilter),taunow); /*/set the tau of the wide-band filter*/
/*/ normalize the gain of the wideband pass filter as 0dB at CF */
dtmp = taunow*TWOPI*(bm->wbfilter.F_shift-bm->bmfilter.F_shift);
wb_gain = pow((1+dtmp*dtmp), bm->wbfilter.Order/2.0);
bm->wbfilter.gain=wb_gain;
ohcout = bm->ohc.run(&(bm->ohc),wbout); /*/ pass the control signal through OHC model*/
/*/ pass the control signal through nonliearity after OHC */
bm->tau = bm->afterohc.run(&(bm->afterohc),ohcout);
/*/ set the tau of the tuning filter */
bm->bmfilter.settau(&(bm->bmfilter),bm->tau);
/*/ Gain Control of the tuning filter */
out1 = pow((bm->tau/bm->TauMax),bm->bmfilter.Order)*x1;
break;
};
out[i] = out1;
};
bm->gfagain.run2(&(bm->gfagain),out,out,length);
return;
}
/* ********************************************************************************************
* * Get TauMax, TauMin for the tuning filter. The TauMax is determined by the bandwidth/Q10
of the tuning filter at low level. The TauMin is determined by the gain change between high
and low level
Also the calculation is different for different species
*/
double Get_tau(int species, double cf,int order, double* _taumax,double* _taumin,double* _taurange)
{
double taumin,taumax,taurange;
double ss0,cc0,ss1,cc1;
double Q10,bw,gain,ratio;
double xcf, x1000;
gain = 20+42*log10(cf/1e3); /*/ estimate compression gain of the filter */
if(gain>70) gain = 70;
if(gain<15) gain = 15;
ratio = pow(10,(-gain/(20.0*order))); /*/ ratio of TauMin/TauMax according to the gain, order */
/*/ Calculate the TauMax according to different species */
switch(species)
{
case 1:
/* Cat parameters: Tau0 vs. CF (from Carney & Yin 88) - only good for low CFs*/
/* Note: Tau0 is gammatone time constant for 'high' levels (80 dB rms) */
/* parameters for cat Tau0 vs. CF */
xcf = cochlea_f2x(species,cf); /* position of cf unit; from Liberman's map */
x1000 = cochlea_f2x(species,1000); /* position for 1 Khz */
ss0 = 6.; cc0 = 1.1; ss1 = 2.2; cc1 = 1.1;
taumin = ( cc0 * exp( -xcf / ss0) + cc1 * exp( -xcf /ss1) ) * 1e-3; /* in sec */
taurange = taumin * xcf/x1000;
taumax = taumin+taurange;
break;
case 0:
/* Human Parameters: From Mike Heinz: */
/* Bandwidths now are based on Glasberg and Moore's (1990) ERB=f(CF,level) equations */
taumax = 1./(TWOPI*1.019*erbGM(cf));
break;
case 9:
/* Universal species from data fitting : From Xuedong Zhang,Ian (JASA 2001) */
/* the Q10 determine the taumax(bandwidths at low level) Based on Cat*/
Q10 = pow(10,0.4708*log10(cf/1e3)+0.4664);
bw = cf/Q10;
taumax = 2.0/(TWOPI*bw);
}; /*/end of switch */
taumin = taumax*ratio;
taurange = taumax-taumin;
*_taumin = taumin;
*_taumax = taumax;
*_taurange = taurange;
return 0;
}
/*/ --------------------------------------------------------------------------------
** Calculate the location on Basilar Membrane from best frequency
*/
double cochlea_f2x(int species,double f)
{
double x;
switch(species)
{
case 0: /*/human */
x=(1.0/0.06)*log10((f/165.4)+0.88);
break;
default:
case 1: /*/cat */
x = 11.9 * log10(0.80 + f / 456.0);
break;
};
return(x);
}
/* --------------------------------------------------------------------------------
** Calculate the best frequency from the location on basilar membrane
*/
double cochlea_x2f(int species,double x)
{
double f;
switch(species)
{
case 0: /*/human
// if((x>35)||(x<0)) error("BM distance out of human range, [in cochlea_x2f(...)]"); */
f=165.4*(pow(10,(0.06*x))-0.88);
break;
default:
case 1: /*/cat */
f = 456.0*(pow(10,x/11.9)-0.80);
break;
};
return(f);
}
/*/ --------------------------------------------------------------------------------
** Calculate the erb at 65dB */
double erb51_65(double CF)
{
double erb1000,erbCf,erb;
erb1000=24.7*(4.37*1000/1000+1);
erbCf=24.7*(4.37*CF/1000+1);
erb=(0.5*erbCf)/(1-(0.38/4000)*erb1000*(65+10*log10(erbCf)-51))+.5*erbCf;
return(erb);
}
/*/ --------------------------------------------------------------------------------
** Calculate the erb using GM's method */
double erbGM(double CF)
{
double erbCf;
erbCf=(1/1.2)*24.7*(4.37*CF/1000+1);
return(erbCf);
}
/*/ ---------------------------------------------------------------------------
** Calculate the delay(basilar membrane, synapse for cat*/
double delay_cat(double cf)
{
/* DELAY THE WAVEFORM (delay buf1, tauf, ihc for display purposes) */
/* Note: Latency vs. CF for click responses is available for Cat only (not human) */
/* Use original fit for Tl (latency vs. CF in msec) from Carney & Yin '88
and then correct by .75 cycles to go from PEAK delay to ONSET delay */
double A0 = 8.13; /* from Carney and Yin '88 */
double A1 = 6.49;
double x = cochlea_f2x(1,cf); /*/cat mapping */
double delay = A0 * exp( -x/A1 ) * 1e-3 - 1.0/cf;
return(delay);
}