annotate wdf/wdf_tl.c @ 0:5242703e91d3 tip

Initial checkin for AIM92 aimR8.2 (last updated May 1997).

/*
    Copyright (c) Applied Psychology Unit, Medical Research Council. 1988, 1989
    ===========================================================================

    Permission to use, copy, modify, and distribute this software without fee 
    is hereby granted for research purposes, provided that this copyright 
    notice appears in all copies and in all supporting documentation, and that 
    the software is not redistributed for any fee (except for a nominal shipping 
    charge). Anyone wanting to incorporate all or part of this software in a
    commercial product must obtain a license from the Medical Research Council.

    The MRC makes no representations about the suitability of this 
    software for any purpose.  It is provided "as is" without express or implied 
    warranty.
 
    THE MRC DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING 
    ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL THE
    A.P.U. BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY 
    DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN 
    AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 
    OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/

/*  
    Acknowledgment:
    ==============

    The source code provided in this file was originally developed by 
    Christian Giguere as part of a Ph.D degree at the Department of
    Engineering of the University of Cambridge from April 1990 to
    November 1993. The code was subsequently adapted under a grant
    from the Hearing Research Trust for full compatibility with 
    AIM Release 6.15.

    Christian Giguere 25/03/94

*/

/*
    =============================================================
    wdf_tl.c
    =============================================================

    Wave digital filter (WDF) implementation of the cochlear 
    transmission line (TL) network.

    Author        : Christian Giguere
    First written : 01st April, 1991
    Last edited   : 18th February, 1994

    References:
    (1) C.Giguere and P.C.Woodland (1994). JASA 95(1): 331-342.
    (2) A.Fettweis (1986). Proceedings IEEE 74(2): 270-327.
    =============================================================
*/

/***** includes *****/

#include <math.h>
#include <stdio.h>
#include "stitch.h"
#include "source.h"
#include "calc.h"
#include "calc_tl.h"
#include "bank_tl.h" 
#include "wdf_tl.h"
#include "wdf_ear.h"
     
/***** defines *****/

#if 0
#define   _DEBUG_
#define   _OVERFLOW_
#endif

#if 0
#define   _IMPULSE_            /* bypasses input wave with a delta impulse */
#endif

#if 0			       /* bypasses input wave with a pure tone and dumps  */
#define   _BMCURVES_           /* the rms value of basilar membrane velocity on stdout */
#endif			       /* (see below for other parameters to set) */

#if 0
#define   _EAR_                /* dumps stapes velocity and eardrum pressure on stdout */
#endif

#define   NUMBER_OF_STATES     (           14 )           /* per TLF segment */
#define   NUMBER_OF_COEFFS     (            0 )           /* per TLF segment */


/***** functions *****/

/******************************** WDF set-up function *************************************
* name:                 function:
*
* WDFilterState()       Set-up function for the WDF implementation of the cochlear network 
*                       (Giguere and Woodland, 1994, Figs. 3 and 10). It generates the 
*                       multiplier coefficients and initializes the states of the filter.
*                       It is called by function ``TLF_GenBank()'' in file ``bank_tl.c'' 
*                       once for each segment of the transmission line starting from the most 
*                       apical segment. It returns a pointer to a structure containing all 
*                       the relevant information for the computation of the current filter 
*                       segment.
*
*                       Note: The code that sets up the WDF implementation of the outer ear 
*                             and middle ear is located in file ``wdf_ear.c''.
******************************************************************************************/

/************************************ WDFilterState() ************************************/

WDFilterState *WDFilter( samplerate, center_frequency, scale_vel, scale_disp, rho, scala_area, 
			 width, qfactor, mass_per_area, seglength, Lt, warping, active, zov, 
			 OHC_gain, OHC_sat )

double  samplerate, center_frequency, scale_vel, scale_disp ;
double  rho, scala_area, width, qfactor, mass_per_area, seglength, Lt ;
int     warping, active ;
double  *zov, OHC_gain, OHC_sat ;
{
    static int  apex = 1 ;
    DeclareNew( WDFilterState *, filter_state ) ;
    CoeffType *cf ;
    double  an, bn, mn, d, fs, fn ; 
    double  Lsn, Ln, Cn, Rn ;
    double  r1, r2, r3, r4, g1, g2, g3, zn, zcn ;

   /***** cochlear parameters for current BM segment *****/
    an = scala_area ;                       /* cross-sectional area of each scala (cm2) */
    bn = width ;                            /* BM width (cm) */
    mn = mass_per_area ;                    /* transversal mass-per-area-of-BM (g/cm2) */
    d = seglength ;                         /* BM segment length (cm) */ 

   /***** frequency warping compensation due to bilinear transform *****/
    fs = samplerate ;
    if ( warping == 0 )
	 fn = center_frequency ;
    else
	 fn = fs / Pi * tan( Pi * center_frequency / fs ) ;

   /***** electrical circuit elements (CGS units) for current segment *****/
    Ln = mn / ( bn * d ) ;                                      /* shunt inductor */
    Cn = 1. / ( ( TwoPi * TwoPi * fn * fn ) * Ln ) ;            /* shunt capacitor */
    Rn = sqrt( Ln / Cn ) / qfactor ;                            /* shunt resistor */
    Lsn =  ( 2. * rho * d ) / an ;                              /* series inductor */


   /***** WDF multiplier coefficients for current TL segment *****/

   /*** initialise filter coeffs ***/
    filter_state->coeffs = ( char * ) NewZeroedArray( CoeffType, NUMBER_OF_COEFFS, 
						      "wdf_tl.c for filter coeffs\n" ) ;
    cf = ( CoeffType * ) filter_state->coeffs ;
    
    /* Adaptor 25 */
    r1 = zcn = 1. / (2. * fs * Cn ) ;                           /* Zcn */
    r2 = Rn ;                                                   /* Zrn */
    r3 = 2. * fs * Ln ;                                         /* Zln */
    r4 = zn = r1 + r2 + r3 ;
    filter_state->gamma251 = r1 / r4 ;
    filter_state->gamma252 = r2 / r4 ;

    /* Adaptor 24 */
    g1 = 1. / zn ;
    if( apex ) {
	g2 = 1. / ( 2. * fs * Lt ) ;                            /* (1 / Zlt) */
	apex = 0 ; 
    }
    else
	g2 = 1. / *zov ;                                        /* Zov (at the base) */     
    g3 = g1 + g2 ;
    filter_state->gamma241 = g1 / g3 ;

    /* Adaptor 23 */
    r1 = 2. * fs * Lsn ;                                        /* Zlsn */
    r2 = 1. / g3 ;
    r3 = *zov = r1 + r2 ;                                       /* Zov */
    filter_state->gamma231 = r1 / r3 ; 

#ifdef _DEBUG_
    fprintf( stderr, "gamma251=%f gamma252=%f\n", filter_state->gamma251
						, filter_state->gamma252 ) ;
    fprintf( stderr, "gamma241=%f            \n", filter_state->gamma241 ) ; 
    fprintf( stderr, "gamma231=%f            \n", filter_state->gamma231 ) ;
#endif

   /***** scaling to convert output to transversal motion of BM segment *****/ 
    filter_state->out_scale_disp = scale_disp * Cn / ( 2. * bn * d ) ;
    filter_state->out_scale_vel = scale_vel / ( 2. * zcn * bn * d ) ;
   
   /***** multiplier coefficients for OHC source *****/
    filter_state->OHC_sat = OHC_sat * scale_disp / filter_state->out_scale_disp ;
    filter_state->OHC_gain = - OHC_gain * Rn * filter_state->OHC_sat / ( 2. * zcn ) ;

#ifdef _BMCURVES_
   /***** initialise rms detection *****/
    filter_state->rms = 0.0 ;
    filter_state->sample = 0 ;
#endif

   /***** initialise filter states *****/
    filter_state->states = ( char * ) NewZeroedArray( StateType, NUMBER_OF_STATES, 
						      "wdf_tl.c for filter states\n" ) ;
    
   /***** is channel active for display ? *****/
    filter_state->active = active ;

   /***** return *****/
    return( filter_state ) ;
}


/*********************************** WDF procedures ***************************************
* name:                 function:
*
* DoWDFdataArray()      WDF realization of the outer ear, middle ear and cochlear networks. 
*                       It is called by function ``tlf_bank_callback()'' in file ``bank_tl.c'' 
*                       once for a specified number of input points. It computes the BM 
*                       displacement or velocity for each segment and keeps track of the 
*                       filter states. 
*
* CloseWDF()            Deletes all private data structures and arrays of the cochlear 
*                       transmission line filter upon completion of filtering. It is called 
*                       by function ``tlf_bank_callbank()'' in file ``bank_tl.c''.
******************************************************************************************/

/********************************** DoWDFdataArray() **********************************/

void DoWDFdataArray( bankInfo, filter_states, input, output, points, channels,
		     ws, eartube_states, ms, tube_segments )
register TLF_BankInfo     *bankInfo ;
register WDFilterState    **filter_states ;
register DataType         *input ;
register DataType         *output ;
register int              points, channels ;
register WaveWDFstate       *ws ;
register EartubeWDFstate    **eartube_states ;
register EarmiddleWDFstate  *ms ;
register int              tube_segments ;
{
    register WDFilterState  *fs ;
    register EartubeWDFstate *ts ;
    register StateType   *st ;   
    register StateType   dn, in, out ; 
    register StateType   a0, an_1, b1, b2, b3, b63, b161, b202, b212 ;
    static   StateType   bn = 0. ;
    register CoeffType   *cf ;
    register int         decimateCount = bankInfo->decimateCount ;
    register int         decimate_factor = bankInfo->decimate_factor ;
    register int         channelActive, pointActive, motion = bankInfo->output_var ;
    register int         channel = -1, tube_segment = -1 ;
#ifdef _BMCURVES_
    static   long        sample = -1 ;
	     int         period = 48 ;                        /* in samples */
	     double      samplerate = 48000. ;                /* samples per sec */
	     double      rise_time = 2.0 ;                    /* sec */
	     double      start_time = 4.0 ;                   /* sec */
#endif
#ifdef _IMPULSE_
    static   int         impulse = 10000 ;
#endif

    /***** update stapedial muscle state *****/

     ( void ) ms->update_proc( ms ) ; 

    /***** for all points ******/

#ifdef _DEBUG_     
     fprintf( stderr, "DoWDFdataArray() = %d points X %d channels\n", points, channels ) ;
#endif     

     while( points-- > 0 ) {

	if( decimateCount == decimate_factor )
	   decimateCount = 0 ;
	pointActive = !decimateCount ;
	decimateCount++ ;

	output += bankInfo->output_chans * pointActive ; 

       /*** for all TL channels in BACKWARD direction for current input point ***/

	while( ++channel < channels ) {
  
	  /* get states and coefficients */
	  fs = *filter_states++ ;   
	  st = ( StateType * ) fs->states ;
	  cf = ( CoeffType * ) fs->coeffs ; 

	  /* Adaptor 25 */
	  st[10] = st[1] ;                                    /* a251 */
	  st[3] = st[2] - st[4] - st[1] ;                     /* b254 = a221 */

	  /* Adaptor 24 */ 
	  st[9] = bn ;                                        /* -a242 */
	  st[5] = fs->gamma241 * ( st[3] + st[9] ) ;          /* b240 */
	  st[6] = st[5] - st[9] ;                             /* b243 = a232 */

	  /* Adaptor 23 */
	  st[7] = st[8] - st[6] ;                             /* b233 */
	  bn = st[7] ;
	}

     /*** input  point ***/

	in = ( StateType ) *input++ ;  

#ifdef _IMPULSE_
	if( impulse != 0 )  {
	   in = ( StateType ) impulse ;
	   impulse = 0 ;
	}
	else 
	   in = ( StateType ) 0 ;
#endif

#ifdef _BMCURVES_
	if( ++sample < ( long ) ( samplerate * rise_time ) )
	    in = 0.28284 * ( double ) sample / ( samplerate * rise_time ) * 
		 sin( TwoPi * ( double ) ( sample % period ) / ( double ) period ) ;
	else
	    in = 0.28284 * sin( TwoPi * ( double ) ( sample % period ) / ( double ) period ) ;
#endif

     /*** computation of outer and middle ear in FORWARD direction for current input point ***/

     /* freefield - external ear */

	/* get states */
	st = ( StateType * ) ws->states ;      

	/* Adaptor 0 */
	st[4] = ws->gamma01 * ( st[3] - in ) ;                /* b00  */
	st[5] = st[4] + in + in ;                             /* b03 = -a11 */ 

	/* Adaptor 2 */
	st[1] = -ws->gamma21 * st[0] ;                        /* b20 */
	st[2] = st[1] ;                                       /* b23 = a12 */

	/* Adaptor 1 */
	st[6] = st[5] - st[2] ;                               /* b13 = a32 */
	an_1 = st[6] ;

     /* concha and auditory canal */

	while( ++tube_segment < tube_segments ) {

	   /* get states */
	   ts = *eartube_states++ ;  
	   st = ( StateType * ) ts->states ;
 
	   /* Adaptor 3 */
	   st[3] = st[4] - an_1 ;                             /* b33 = a41 */

	   /* Adaptor 5 */
	   st[0] = -ts->gamma51 * st[1] ;                     /* b50 */
	   st[2] = st[0] + st[1] ;                            /* b53 = a42 */
 
	   /* Adaptor 4 */
	   st[5] = st[3] - st[2] ;                            /* -a42+a41 */
	   st[6] = ts->gamma41 * st[5] ;                      /* b40 */
	   st[7] = st[6] + st[2] ;                            /* b43 = a61 */

	   /* Adaptor 6 */
	   st[8] = st[9] - st[7] ;                            /* b63 = a32 */
	   an_1 = b63 = st[8] ;
	}

     /* middle ear */

	/* get states */
	st = ( StateType * ) ms->states ;      

	/* Adaptor 9 */
	st[31] = st[30] - st[29] ;                           /* b94 = a81 */

	/* Adaptor 8 */
	st[32] = st[35] - st[31] ;                           /* a83-a81 */
	st[33] = -ms->gamma81 * st[32] - ms->gamma82 * st[35] ;   /* b80 */ 
	st[34] = st[33] + st[35] ;                           /* b84 = a71 */

	/* Adaptor 7 */
	st[36] = -st[34] - an_1 ;                            /* b73 = a102 */

	/* Adaptor 13 */
	st[15] = -st[14] ;                                   /* b133 = a122 */

	/* Adaptor 12 */
	st[16] = st[15] + st[18] ;                           /* a122-a121 */
	st[17] = -ms->gamma121 * st[16] ;                    /* b120 */
	st[19] = st[17] + st[15] ;                           /* b123 = a112 */

	/* Adaptor 14 */
	st[21] = -ms->gamma141 * st[20] ;                    /* b140 */
	st[22] = st[21] + st[20] ;                           /* b143 = a111 */

	/* Adaptor 15 */
	st[23] = -st[24] ;                                   /* b153 = a113 */

	/* Adaptor 11 */
	st[25] = -st[22] - st[19] - st[23] ;                 /* b114 = a101 */

	/* Adaptor 10 */
	st[26] = st[25] - st[36] ;                           /* -a102+a101 */
	st[27] = ms->gamma101 * st[26] ;                     /* b100 */
	st[28] = st[27] + st[36] ;                           /* b103 = a161 */

	/* Adaptor 17 */
	st[10] = st[12] - st[11] ;                           /* b174 = a162 */
	
	/* Adaptor 16 */
	st[13] = -st[28] - st[10] ;                          /* b163 = a181 */

	/* Adaptor 19 */
	st[5] = -st[6] ;                                     /* b193 = a182 */

	/* Adaptor 18 */
	st[7] = st[13] - st[5] ;                             /* -a182+a181 */
	st[8] = ms->gamma181 * st[7] ;                       /* b180 */
	st[9] = st[8] + st[5] ;                              /* b183 = a203 */

	/* Adaptor 20 */
	st[40] =  - bn / ms->ratio ;                         /* a202 */
	st[3] = -st[40] - st[9] ;                            /* b204 = a212 */

	/* Adaptor 21 */
	st[0] = st[2] - st[1] - st[3] ;                      /* b214 = a222 */


     /*** computation of outer and middle ear in BACKWARD direction for current input point ***/

     /* middle ear */

	/* Adaptor 22 */
	a0 = st[4] + st[0] ;                                  /* a220 */
	st[4] = st[4] - ms->gamma221 * a0 ;                   /* b221 */
	b2 = -st[4] - a0 ;                                    /* b222 = a214 */

	/* Adaptor 21 */
	a0 = b2 - st[0] ;                                     /* a210 */
	st[1] = st[1] - ms->gamma211 * a0 ;                   /* b211 */
	b212 = st[3] - ms->gamma212 * a0 ;                    /* b212 = a204 */
	st[2] = -st[1] - b212 - b2 ;                          /* b213 */

	/* Adaptor 20 */
	a0 = b212 - st[3] ;                                   /* a200 */
	b202 = st[40] - ms->gamma202 * a0 ;                   /* b202 */
	b3 = ms->gamma201 * a0 - b202 - b212 ;                /* b203 = a183 */

	/* Adaptor 18 */
	b2 = st[8] + b3 ;                                     /* b182 = a193 */
	b1 = b2 - st[7] ;                                     /* b181 = a163 */

	/* Adaptor 19 */
	st[6] = st[6] + ms->gamma191 * ( st[5] - b2 ) ;       /* b191 */

	/* Adaptor 16 */
	b161 = st[28] + ms->gamma161 * ( st[13] - b1 ) ;      /* b161 = a103 */
	b2 = -b161 - b1 ;                                     /* b162 = a174 */

	/* Adaptor 17 */
	a0 = b2 - st[10] ;                                    /* a170 */
	st[11] = st[11] - ms->gamma171 * a0 ;                 /* b171 */
	st[12] = ms->gamma172 * a0 - b2 - st[11] ;            /* b173 */
 
	/* Adaptor 10 */
	st[37] = st[27] + b161 ;                              /* b102 = a73 */
	st[38] = st[37] - st[26] ;                            /* b101 = a114 */

	/* Adaptor 11 */
	a0 = st[38] - st[25] ;                                /* a110 */
	b1 = st[22] - ms->gamma111 * a0 ;                     /* b111 = a143 */
	b2 = st[19] - ms->gamma112 * a0 ;                     /* b112 = a123 */
	b3 = -b1 -b2 - st[38] ;                               /* b113 = a153 */

	/* Adaptor 15 */
	st[24] = st[24] + ms->gamma151 * ( st[23] - b3 ) ;    /* b151 */

	/* Adaptor 14 */
	st[20] = st[21] + b1 ;                                /* b142 */

	/* Adaptor 12 */
	b2 = st[17] + b2 ;                                    /* b122 = a133 */
	st[18] = b2 + st[16] ;                                /* b121 */

	/* Adaptor 13 */
	st[14] = st[14] + ms->gamma131 * ( st[15] - b2 ) ;    /* b131 */

	/* Adaptor 7 */
	b1 = st[34] + ms->gamma71 * ( st[36] - st[37] ) ;     /* b71 = a84 */
	st[39] = -b1 - st[37] ;                               /* b72 = a63 */

	/* Adaptor 8 */
	st[35] = st[33] + b1 ;                                /* b83 */
	b1 = st[35] + st[32] ;                                /* b81 = a94 */

	/* Adaptor 9 */
	a0 = b1 - st[31] ;                                    /* a90 */
	st[29] = st[29] - ms->gamma91 * a0 ;                  /* b91 */
	st[30] = -st[29] + ms->gamma92 * a0 - b1 ;            /* b93 */

     bn = st[39] ;

     /* concha and auditory canal */

	while( tube_segment-- > 0 ) {

	   /* get states */
	   ts = *--eartube_states ;   
	   st = ( StateType * ) ts->states ;

	   /* Adaptor 6 */
	   b1 = st[7] - ts->gamma61 * ( bn - st[8] ) ;        /* b61 = a43 */
	   st[9] = -b1 - bn ;                                 /* b62 */

	   /* Adaptor 4 */
	   b2 = st[6] + b1 ;                                  /* b42 = a53 */
	   b1 = b2 - st[5] ;                                  /* b41 = a33 */

	   /* Adaptor 5 */
	   st[1] = st[0] + b2 ;                               /* b52 */

	   /* Adaptor 3 */
	   st[4] = ts->gamma31 * ( st[3] - b1 ) - st[4] ;     /* b31 */
	   bn = -st[4] - b1 ;                                 /* b32 = a63 */
	}

     /* freefield - external ear */

	/* get states */
	st = ( StateType * ) ws->states ;      

	/* Adaptor 1 */
	b1 = ws->gamma11 * ( st[6] - bn ) - st[5] ;           /* b11 = -a03 */
	b2 = -b1 - bn ;                                       /* b12 = a23 */

	/* Adaptor 2 */
	st[0] = st[1] + b2 + st[0] ;                          /* b61 */

	/* Adaptor 0 */
	st[3] = b1 - st[4] + st[3] ;                          /* -b01 + p */

     an_1 = - b202 * ms->ratio ;    


	/* for all channels in FORWARD direction for current input point */

	 while( channel-- > 0 ) {

	   /* get states and coefficients */
	   fs = *--filter_states ;   
	   st = ( StateType * ) fs->states ;
	   cf = ( CoeffType * ) fs->coeffs ;
	   channelActive = fs->active ;

	   /* Adaptor 23 */
	   st[8] = fs->gamma231 * ( st[7] - an_1 ) - st[8] ;    /* b231 */
	   b2 = an_1 + st[8] ;                                  /* -b232 = -a243 */

	   /* Adaptor 24 */
	   an_1 = b2 - st[5] ;                                  /* -b242 */
	   b1 = -an_1 - st[3] - st[9] ;                         /* b241 = a254 */
	  
	   /* Adaptor 25 */
	   a0 = b1 - st[3] ;                                    /* a250 */
	   st[1] = st[10] - fs->gamma251 * a0 ;                 /* b251 */
	   st[2] = fs->gamma252 * a0 - st[4] - st[1] - b1 ;     /* b253 */
	   dn = st[1] + st[10] ; 
	   in = st[1] - st[10] ;

	   /* output storage */
	   if( channelActive * pointActive ) {

	      if( motion == DISPLACEMENT ) 
		 out = fs->out_scale_disp * dn ;
	      else 
		 out = fs->out_scale_vel * in ;
#ifdef _OVERFLOW_
	      if( out > _MaxOutput_ || out < _MinOutput_ )
		 fprintf( stderr, "Overflow error in BMM output=%e\n", ( double ) out ) ;     
#endif
	      *(--output) = ( DataType ) ( out ) ; 
	   }

	   /* OHC nonlinear voltage source */
	      if( dn < 0. )
		 dn = -dn ;
	      st[4] = fs->OHC_gain * in / ( fs->OHC_sat + dn ) ;  /* -Vn(ohc) */

#ifdef _BMCURVES_
	   if( sample > ( long ) ( start_time * samplerate ) ) { 
	      in = fs->out_scale_vel * in ;                     
	      fs->rms = fs->rms + in*in ; 
	      fs->sample++ ;
	   }
#endif

	}
	bn = -an_1 ;                                     /* apical boundary condition */
	output += bankInfo->output_chans * pointActive ;

#ifdef _EAR_
	st = ( StateType * ) ms->states ;
	printf( "stapes=%e  eardrum=%e\n", 0.5 * ( st[40] - b202 ), 0.5 * ( st[39] + b63 ) ) ; 
#endif

    }

    return ;
}


/******************************* CloseWDF() *************************************/

void CloseWDF( states, channels, bank )
  register WDFilterState      **states ;
  register int                channels ;
  register TLF_BankInfo       *bank ;
{
  int channel ;

  for( channel = 0 ; channel < channels ; channel++ ) {
#ifdef _BMCURVES_
     printf( "%d  %e\n", channel+1, sqrt(states[channel]->rms / states[channel]->sample) ) ;
#endif
     Delete( states[channel]->states ) ;
     Delete( states[channel]->coeffs ) ;
  }
  Delete( states ) ;
  
  Delete( bank ) ;

  return ;

}