anmodel.cpp

/** \file anmodel.cpp
 * \brief This is the main file for the AN model.
 * <p>the output is stored in filterout, synapseout
 * <p>The stimulus is read from the file(default filename is "click")
 * <p>The program initializes the filter bank model first(according to the parameters
 * specified in the commandline).
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <iostream.h>
#include <fstream.h>

#include "cmpa.hpp"
#include "complex.hpp"

#define NSTIMMAX 100000 // enough points for durations up to 250 msec
#define NUMCHANS 64 // Number of fibers in the population to be simulated
#define MAXSPIKES 10000
double stim[NSTIMMAX];
double sout[NSTIMMAX];
double sptime[MAXSPIKES]; //spike time storage

/// This structure defines the parameters needed for the model
typedef struct {
  char wavefile[100]; //filename for stimulus waveform
  int banks; //number of filters in filterbank
  double delx;  //"delta x", or spacing of "IHC"s along BM (linearly spaced IHCs in terms of location along BM
				//  results in approximately logarithmic spacing in terms of frequency.
  int stim;
  double tdres; //time step (time domain resolution)
  int nstim,nrep; // nstim=number of pts in stimulus; nrep is number of reps of stimulus
  int runmode;
  double freqt,spl,dur,reptim;    //for tone and other stim
                                  //noise(2),click(3)
  double F1,F2,spl1,spl2;         //for two tone suppression(8)
  double *buf;
  double rf;                      //rise & fall
  int rfmethod;// rise & fall method
  double cfhi,cflo,cfinc;
  double splhi,spllo,splinc;
  double freqthi,freqtlo,freqtinc;
} T_stim;

int parsecommandline(TAuditoryNerve *pan,T_stim *ptm,int argc, char* argv[]);
int synapse(TAuditoryNerve *pan,int nstim,int nrep,double *sout,double *sptime);

int main(int argc, char* argv[])
{
  static T_stim tm;                    // structure to store the input parameters
  static TAuditoryNerve an[NUMCHANS];              // filter bank(Maximum 64)
  double xlo,xhi,delx,xcenter;

  int species = 1; //CAT
  int icenter;
  int icf;
  int nspikes;

  double x,cf;
  FILE *fpstim,*fpbm,*fpsout;          //open these file when run
  double buf; //tempory value to store the stimulus
  long flag,nowstim,totalstim;

  an[0].tdres = tm.tdres = 10e-6;      //default sample rate = 500kHz
  an[0].cf = 1000;                     //default CF      : 1000 Hz
  an[0].spont = 50;                    //     spont rate : 50
  strcpy(tm.wavefile,"click");         //input file      : click
  tm.reptim = 0.020;                   //repetition time  : 20 msec
  tm.delx = 0.05;                      //distance between two filters along BM : 0.05mm
  tm.cflo = -1;
  tm.cfhi = -1;
  tm.banks = -1;                       //default         : no filters
  tm.nrep = 0;

  parsecommandline(&an[0],&tm,argc,argv); // get the parameters from command line

  if(tm.banks == -1) //this is for a simulation of a single anditory nerve fiber
    {
      tm.banks = 1;
      xcenter = cochlea_f2x(species,an[0].cf); //convert from CF to BM location
      delx = 0;
    }
  else if(tm.cflo == -1) //this is for AN fiber population
    {
      xcenter = cochlea_f2x(species,an[0].cf);
      delx = tm.delx;
    }
  else
    {
      if(tm.cfhi<tm.cflo) error("cfhi should be greater than cflo");
      xlo = cochlea_f2x(species,tm.cflo);
      xhi = cochlea_f2x(species,tm.cfhi);
      xcenter = (xhi -xlo)/2. + xlo;
      delx = (xhi-xlo)/(tm.banks-1);
      an[0].cf = cochlea_x2f(species,xcenter);
    };

  printf("\n Filter Bank: center CF=%f, xcenter=%f(mm from apex), number of fibers= %d\n",an[0].cf,xcenter,tm.banks);

  icenter = tm.banks/2;

  nowstim = 0;
  totalstim = (long)(tm.reptim/an[0].tdres);
  buf = 0.0;

  // ******************now begin to pass the stimulus to the model************************
  // Read in the stimulus
  for(nowstim=0; nowstim<NSTIMMAX; nowstim++) stim[nowstim] = 0.0;
  printf("\nReading in stimulus \n");
  fpstim = fopen(tm.wavefile,"r");
  nowstim=0;
  while((fscanf(fpstim,"%le",&stim[nowstim]) != EOF) & nowstim<NSTIMMAX)
    nowstim++;
  printf("%d points read in\n",nowstim);
  if(totalstim < nowstim) totalstim=nowstim;
  fclose(fpstim);
  if(nowstim>=NSTIMMAX) error("Entire stim not read in: increase NSTIMMAX");
  if(totalstim>=NSTIMMAX) error("Entire reptim too long: increase NSTTMMAX");

  // construct the model
  for(icf = 0; icf < tm.banks; icf++) /* Loop through the filters */
    {
      x = xcenter + (icf - icenter)*delx;
      cf = cochlea_x2f(species,x); // from BM location to CF
      an[icf].cf = cf;
      an[icf].tdres = an[0].tdres;
      an[icf].spont = an[0].spont;
      an[icf].construct();
      an[icf].init(cf,an[icf].spont);
      printf("icf = %d out of %d x= %f cf= %f\n",icf,tm.banks,x,cf);
    };

  flag = 0; //if = 1,don't read from stimulus file but set zero as input
  nowstim = 0;
  fpbm = fopen("filterout","w");
  fpsout = fopen("synapseout","w");
  while(nowstim<totalstim)
    {
	  buf = stim[nowstim];
      fprintf(fpbm,"%g",nowstim*an[0].tdres);
      fprintf(fpsout,"%g",nowstim*an[0].tdres);
      for(icf = 0; icf < tm.banks; icf++) // Loop through the filters
		{
		  an[icf].run(buf);
		  fprintf(fpbm," %g",an[icf].bmbuf);  // Write filter response out to file
		  fprintf(fpsout," %g",an[icf].sbuf);  // Write synapse output to file
	  };
      fprintf(fpbm,"\n");
      fprintf(fpsout,"\n");
	  sout[nowstim] = an[0].sbuf;
      nowstim++;
    }; //end of while
  fclose(fpbm);
  fclose(fpsout);
  printf("\n AN filter bank simulation is complete. \n");
  if((tm.nrep>0)&&(tm.banks==1))
	{
	  printf("\n Generate spike times and the PST histogram(write out to files pstplot and sptime...\n");
      nspikes = synapse(&an[0],totalstim,tm.nrep,sout,sptime);
      double bin = 0.0001;
	  int nbins = (int)(tm.reptim / bin + 1);
	  int *pst = new int[nbins];
	  int i;
	  /* Construct PST HISTOGRAM, Interval hist, Period hist.  */
	  for(i = 0; i < nbins; i++) pst[i] = 0;
	  for(i = 0; i < nspikes; i++)
      {
      	/* both sptime and bin are in secs */
	    int ipst = (int)(fmod(sptime[i],tm.tdres*totalstim) / bin);
        pst[ipst] = pst[ipst] + 1;
      }
 	  fstream *fp = new fstream("pstplot",ios::out);
      for(i = 0; i < nbins; i++)
	  {
 		  	(*fp)<<(i*bin)<<' '<<pst[i]<<"\n";
      };
      delete fp;
      delete pst;

 	  fp = new fstream("sptime",ios::out);
      for(i = 0; i < nspikes; i++)
	  {
 		  	(*fp)<<i<<' '<<sptime[i]<<"\n";
      };
      delete fp;
	};
};

// ----------------------------------------------------------------
/**
  * This function reads the command line options and store them in the structure
*/
int parsecommandline(TAuditoryNerve *pan,T_stim *ptm,int argc, char* argv[])
{
  int i;
  int needhelp = 0;
  char *para;

  i = 1;
  if(argc ==1 ) needhelp = 1;
  while(i<argc)
  { para = argv[i];
    if((*para)=='-')
      { para++;
	if(strcmp(para,"tdres")==0)
	  {
	    pan->tdres = double(atof(argv[++i]));
	    ptm->tdres = pan->tdres;
	  }
	else if(strcmp(para,"cf")==0)
	  pan->cf = double(atof(argv[++i]));
	else if(strcmp(para,"spont")==0)
	  pan->spont = double(atof(argv[++i]));

	else if(strcmp(para,"fibers")==0)
	  ptm->banks = atoi(argv[++i]);
	else if(strcmp(para,"reptim")==0)
	  ptm->reptim = double(atof(argv[++i]));
	else if(strcmp(para,"delx")==0)
	  ptm->delx = double(atof(argv[++i]));
	else if(strcmp(para,"trials")==0)
	  ptm->nrep = atoi(argv[++i]);
	else if(strcmp(para,"cfhi")==0)
	  ptm->cfhi = double(atof(argv[++i]));
	else if(strcmp(para,"cflo")==0)
	  ptm->cflo = double(atof(argv[++i]));
	else if(strcmp(para,"wavefile")==0)
	  {
	    ptm->stim = 11;
	    strcpy(ptm->wavefile,argv[++i]);
	  }
	else if(strcmp(para,"help")==0)
	  needhelp = 1;
	else
	  { printf("\nUnkown parameters --> %s",para); needhelp = 1; break; };
      }
    else { printf("\nUnkown parameters --> %s",para); needhelp = 1; break; };
    i++;
  };
  if(needhelp==1)
    {
      printf("\n This program accepts following parameters:\n");

      printf("\n -cf #(1000)   --> character frequency of the AN model fiber (or center CF for filter bank)");
      printf("\n -spont #(50)  --> spontaneous rate of the fiber");
      printf("\n -tdres #(10e-6)--> time domain resolution(in seconds)");

      printf("\n -fibers #(1) --># of filters, use this option with cf,cflo,cfhi,delx");
      printf("\n -wavefile filename(click) --> specify the stimulus wavefile name(click)");
      printf("\n -delx #(0.05mm) --> distance between filters along basilar membrane");
      printf("\n -reptim #(0.02) --> Total duration of stimulation(20msec)");
      printf("\n -cfhi #(-1)   --> highest cf in population (specify cfhi,cflo will recalculate cf&delx");
      printf("\n -cflo #(-1)   --> lowest cf in population");
      printf("\n -trials #(0)  --> spike generation trials (this is valid only if # of fibers=1)");
      printf("\n");
      exit(1);
    };
return(0);
};

// --------------------------------------------------------------------------------
/** this routine generates spike times according to *sout with nrep representation
   put the spikes time into buffer *spikes, nspikes
*/
int synapse(TAuditoryNerve *pan,int nstim,int nrep,double *sout,double *sptime)
{
  int i,j,isp;
  int nspikes = 0;
  /* sout[] contains prob of spike vs. time - run through it nrep times to look for spikes */
  isp = 0;
  sptime[0] = 0.; /* start out with a spike at t=0 on 1st rep */
  pan->sg->init(sout[0]); //initialize the spike generator
  for(i = 0; i < nrep; i++)
    {
      pan->sg->init(sout[0]);
      for( j = 0; j < nstim; j++)
	{
	  if(pan->sg->run(sout[j])==1)
	    {
	      sptime[isp] = pan->sg->Get_rtime();  /* leave sptime in sec */
	      isp++;
		  if(isp>MAXSPIKES) error("\nToo much spikes generated: increase MAXSPIKES\n");
	    }
	}
      nspikes = isp;
    };
  return nspikes;
};