twocadap.c

Auditory Simulation Development Computing System version 1.5.2, is based upon a unified re-interp

/*******
 *
 * TwoCAdapt.c
 *
 * This is a simple test harness to make sure that things are working.
 * This harness calculates the two component adaptation functions for the hair
 * cell model.
 * It uses an approximation which can only be used with smooth spike
 * probabilities i.e. Meddis 86 IHC model.
 *
 * by Lowel O'Mard 17-2-93.
 *
 ********/

#include <stdlib.h> 
#include <stdio.h>
#include <math.h>
#include <time.h>

#include "CRLHeaders.h"

/******************************************************************************/ 
/****************************** Constant definitions **************************/ 
/******************************************************************************/ 

#define	PARAMETERS_FILE		"TwoCAdapt.par" /* Name of paramters file. */
#define	NUM_CHANNELS		1			/* Number of signal channels. */
#define	CHANNEL				0			/* Channel for testing. */

/******************************************************************************/ 
/****************************** Global variables ******************************/ 
/******************************************************************************/

char	outputFile[MAXLINE], stParFile[MAXLINE], pEParFile[MAXLINE];
char	bMParFile[MAXLINE], rPParFile[MAXLINE], hCParFile[MAXLINE];
char	sGParFile[MAXLINE];

char	stModuleName[MAXLINE], pEModuleName[MAXLINE], bMModuleName[MAXLINE];
char	rPModuleName[MAXLINE], hCModuleName[MAXLINE], sGModuleName[MAXLINE];

int		numReadings;
double	initialIntensity, intensityIncrement;
double	rampInterval;
double	shortTermT1, shortTermPeriod, rapidAdaptPeriod;
double	binWidth;

/******************************************************************************/ 
/****************************** Functions and subroutines *********************/ 
/******************************************************************************/

/****************************** ReadParsFromFile ******************************/

/*
 * This program reads a specified number of parameters from a file.
 * It expects there to be one parameter per line.
 */
 
void
ReadParsFromFile(char *fileName)
{
	char	line[MAXLINE];
	FILE	*fp;
	
	if ((fp = fopen(fileName, "r")) == NULL) {
		NotifyError("ReadTestPars: Cannot open data file '%s'.\n", fileName);
		exit(1);
	}
	printf("Reading parameters from file: %s\n", fileName);
	Init_ParFile();
	GetPars_ParFile(fp, "%s", outputFile);
	GetPars_ParFile(fp, "%s %s", stParFile, stModuleName);
	GetPars_ParFile(fp, "%s %s", pEParFile, pEModuleName);
	GetPars_ParFile(fp, "%s %s", bMParFile, bMModuleName);
	GetPars_ParFile(fp, "%s %s", rPParFile, rPModuleName);
	GetPars_ParFile(fp, "%s %s", hCParFile, hCModuleName);
	GetPars_ParFile(fp, "%s %s", sGParFile, sGModuleName);
	GetPars_ParFile(fp, "%d", &numReadings);
	GetPars_ParFile(fp, "%lf", &initialIntensity);
	GetPars_ParFile(fp, "%lf", &intensityIncrement);
	GetPars_ParFile(fp, "%lf", &rampInterval);
	GetPars_ParFile(fp, "%lf", &shortTermT1);
	GetPars_ParFile(fp, "%lf", &shortTermPeriod);
	GetPars_ParFile(fp, "%lf", &rapidAdaptPeriod);
	GetPars_ParFile(fp, "%lf", &binWidth);
	fclose(fp);
	Free_ParFile();
	
}

/****************************** ReadParsFromFile ******************************/

/*
 * This subroutine calculates the best fit line parameters for a set of data
 * points, x,y.
 * It returns the parameters as arguments.
 */

void
LinReg(double *y0, double *gradient, ChanData y[], double dx, ChanLen length)
{
	int     i;
	double  a, b, c, e, f, x;

	b = c = e = f = 0.0;
	a = (double) length;
	for (i = 0; i < length; i++) {
		x = dx * i;	
		b += x;
		c += x * x;
		e += y[i];
		f += x * y[i];
	}
	*gradient = (f - b * e / a) / (c - b * b / a);
	*y0 = (e - b * *gradient) / a;

} /* LinReg */

/****************************** Calc2CompAdapt ********************************/

/*
 * This analysis routine calculates the two component adaptation constants.
 * It expects the "work" EarObject to be connected to a post stimulus time
 * histogram.
 * The algorithm expects the signal to have a long length.
 * It assumes that the highest peak is the onset peak.
 * Ass is calculated over the last 10% of the signal.
 * This routine destroys its output signal when it has finished using it.
 *
 */

BOOLEAN
Calc2CompAdapt(double *Tst, double *Tr, EarObjectPtr data, double shortTermT1,
  double shortTermPeriod, double rapidAdaptPeriod)
{
	int		chan;
	double	Ast, Ass[MAX_CHANNELS];
	double	gradient, constant, maxValue, yRT1, yRT2, dt, duration, deltaY;
	register	ChanData	*inPtr, *outPtr;
	ChanLen	sT1Indx, sTPeriodIndx, i, rT1Indx, rT2Indx;
	
	if (data == NULL) {
		NotifyError("Calc2CompAdapt_ExpAnalysis: EarObject not initialised.");
		return(FALSE);
	}
	if (!CheckPars_SignalData(data->inSignal)) {
		NotifyError("Calc2CompAdapt_ExpAnalysis: Input signal not "\
		  "correctly set.");		
		return(FALSE);
	}
	dt = data->inSignal->dt;
	sT1Indx = (ChanLen) (shortTermT1 / dt);
	if (sT1Indx > data->inSignal->length) {
		NotifyError("Calc2CompAdapt_ExpAnalysis: Invalid short term T1 (%g ms).",
		  MSEC(shortTermT1));
		return(FALSE);
	}
	sTPeriodIndx = (ChanLen) (shortTermPeriod / dt);
	if (sTPeriodIndx + sT1Indx > data->inSignal->length) {
		NotifyError("Calc2CompAdapt_ExpAnalysis: Invalid short term period "\
		  "(%g ms).", MSEC(shortTermPeriod));
		return(FALSE);
	}
	SetProcessName_EarObject(data, "Two component adaptation calculation.");
	if (!InitOutSignal_EarObject(data, data->inSignal->numChannels,
	  data->inSignal->length, data->inSignal->dt)) {
		NotifyError("CalcPSTH_ExpAnalysis: Cannot initialise output channels.");
		return(FALSE);
	}	
	duration = GetDuration_SignalData(data->inSignal);
	CalcAverages_GenAnalysis(Ass, data->inSignal, duration * 0.9,
	  duration);
	for (chan = 0; chan < data->inSignal->numChannels; chan++) {
		inPtr = data->inSignal->channel[chan];
		outPtr = data->outSignal->channel[chan];
		for (i = 0; i < data->inSignal->length; i++) {
			deltaY = *inPtr++ - Ass[chan];
			*outPtr++ = log(deltaY * deltaY);
		}
		LinReg(&constant, &gradient, &data->outSignal->channel[chan][sT1Indx],
		  dt, sTPeriodIndx);
		*Tst = -2.0 / gradient;
		Ast = exp(constant / 2.0);
		inPtr = data->inSignal->channel[chan];
		for (i = 0, maxValue = 0.0, inPtr = data->inSignal->channel[chan]; i <
		  data->inSignal->length; i++, inPtr++)
			if (*inPtr > maxValue) {
				maxValue = *inPtr;
				rT1Indx = i;
			}
		rT2Indx = rT1Indx + (ChanLen) (rapidAdaptPeriod / dt);
		yRT1 = data->inSignal->channel[chan][rT1Indx] - Ast * exp(-(rT1Indx *
		  dt) / *Tst);
		yRT2 = data->inSignal->channel[chan][rT2Indx] - Ast * exp(-(rT2Indx *
		  dt) / *Tst);
		*Tr = (rapidAdaptPeriod) / (log(yRT1 - Ass[chan]) - log(yRT2 -
		  Ass[chan]));
	}
    if (data->updateCustomersFlag)
    	UpdateCustomers_EarObject(data);
    return(TRUE);
	
}

/******************************************************************************/ 
/****************************** Main Body *************************************/ 
/******************************************************************************/ 

int main()
{
	int		i;
	double	intensity, Tst, Tr;
	FILE	*fp;
	EarObjectPtr	stimulus = NULL, pEFilter = NULL, bMFilter = NULL;
	EarObjectPtr	recpPotn = NULL, hairCell = NULL, spikeGeneration = NULL;
	EarObjectPtr	pSTHistogram = NULL;
	
	printf("Starting Test Harness...\n");
	
	ReadParsFromFile(PARAMETERS_FILE);
	printf("This test routine calculates the two component adaptation ");
	printf("functions for the Meddis 1 hair cell.\n");
	printf("It uses the approximation in Meddis 1988 (JASA).\n");
	printf("A %s stimulus of frequency ramped with a rise time of %lg\n",
	  stModuleName, MSEC(rampInterval));
	printf("is used.  %d readings are taken, with stimulus values increasing\n",
	  numReadings);
	printf("from %lg dB (SPL) in increments of %lg dB.\n", initialIntensity,
	  intensityIncrement);
	printf("The post stimulus time histogram bin width used for the\n");
	printf("adaptation approximation is %lg ms.\n\n", MSEC(binWidth));
	
	/* Initialise EarObjects */
	
	if ((stimulus = Init_EarObject(stModuleName)) == NULL)
		exit(1);
	if ((pEFilter = Init_EarObject(pEModuleName)) == NULL)
		exit(1);
	if ((bMFilter = Init_EarObject(bMModuleName)) == NULL)
		exit(1);
	if ((recpPotn = Init_EarObject(rPModuleName)) == NULL)
		exit(1);
	if ((hairCell = Init_EarObject(hCModuleName)) == NULL)
		exit(1);
	if ((spikeGeneration = Init_EarObject(sGModuleName)) == NULL)
		exit(1);
	pSTHistogram = Init_EarObject("null");
    
    /* Set up Connections */

	ConnectOutSignalToIn_EarObject( stimulus, pEFilter );
	ConnectOutSignalToIn_EarObject( pEFilter, bMFilter );
	ConnectOutSignalToIn_EarObject( bMFilter, recpPotn );
	ConnectOutSignalToIn_EarObject( recpPotn, hairCell );
	ConnectOutSignalToIn_EarObject( hairCell, spikeGeneration );
 	ConnectOutSignalToIn_EarObject( spikeGeneration, pSTHistogram );
	
	/* Initialise modules. */

	printf("Module parameters...\n\n" );

	if (!DoFun1( ReadPars, stimulus, stParFile))
		exit(1);
	DoFun( PrintPars, stimulus );

	if (!DoFun1( ReadPars, pEFilter, pEParFile))
		exit(1);
	DoFun( PrintPars, pEFilter );

	if (!DoFun1( ReadPars, bMFilter, bMParFile))
		exit(1);
	DoFun( PrintPars, bMFilter );

	if (!DoFun1( ReadPars, recpPotn, rPParFile))
		exit(1);
	DoFun( PrintPars, recpPotn );
	
	if (!DoFun1( ReadPars, hairCell, hCParFile))
		exit(1);
	DoFun( PrintPars, hairCell );

	if (!DoFun1( ReadPars, spikeGeneration, sGParFile))
		exit(1);
	DoFun( PrintPars, spikeGeneration );
			
	/* Start main process. */

	printf("\n");
	if ((fp = fopen(outputFile, "w")) == NULL) {
		fprintf(stderr, "TwoCAdapt: Cannot open file '%s'.\n", outputFile);
		exit(1);
	}
	for (i = 0; i < numReadings; i++) {
		printf("Calculating reading %d of %d...\n", i + 1, numReadings);
		intensity = initialIntensity + (i * intensityIncrement);
		DoFun1(SetIntensity, stimulus, intensity);
		DoProcess(GenerateSignal, stimulus);
		if (!stimulus->outSignal->rampFlag )
			RampUpOutSignal_Ramp(stimulus, Sine_Ramp, rampInterval );
		DoProcess(RunModel, pEFilter);
		DoProcess(RunModel, bMFilter);
		DoProcess(RunModel, recpPotn);
		DoProcess(RunModel, hairCell);
		DoProcess(RunModel, spikeGeneration);
		CalcPSTH_ExpAnalysis(pSTHistogram, binWidth);
		Calc2CompAdapt(&Tst, &Tr, pSTHistogram, shortTermT1, shortTermPeriod,
		  rapidAdaptPeriod);
		fprintf(fp, "%10.5lf \t %10.5lf \t %10.5lf\n", intensity, Tst, Tr);
	}
	fclose(fp);

	FreeAll_EarObject();
	
	printf("Finished test.\n");
	return(0);
	
}