tomwalters@0: /*
tomwalters@0:     Copyright (c) Applied Psychology Unit, Medical Research Council. 1988, 1989
tomwalters@0:     ===========================================================================
tomwalters@0: 
tomwalters@0:     Permission to use, copy, modify, and distribute this software without fee 
tomwalters@0:     is hereby granted for research purposes, provided that this copyright 
tomwalters@0:     notice appears in all copies and in all supporting documentation, and that 
tomwalters@0:     the software is not redistributed for any fee (except for a nominal shipping 
tomwalters@0:     charge). Anyone wanting to incorporate all or part of this software in a
tomwalters@0:     commercial product must obtain a license from the Medical Research Council.
tomwalters@0: 
tomwalters@0:     The MRC makes no representations about the suitability of this 
tomwalters@0:     software for any purpose.  It is provided "as is" without express or implied 
tomwalters@0:     warranty.
tomwalters@0:  
tomwalters@0:     THE MRC DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING 
tomwalters@0:     ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL THE
tomwalters@0:     A.P.U. BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY 
tomwalters@0:     DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN 
tomwalters@0:     AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 
tomwalters@0:     OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: /*  
tomwalters@0:     Acknowledgment:
tomwalters@0:     ==============
tomwalters@0: 
tomwalters@0:     The source code provided in this file was originally developed by 
tomwalters@0:     Christian Giguere as part of a Ph.D degree at the Department of
tomwalters@0:     Engineering of the University of Cambridge from April 1990 to
tomwalters@0:     November 1993. The code was subsequently adapted under a grant
tomwalters@0:     from the Hearing Research Trust for full compatibility with 
tomwalters@0:     AIM Release 6.15.
tomwalters@0: 
tomwalters@0:     Christian Giguere 25/03/94
tomwalters@0: 
tomwalters@0: */
tomwalters@0: 
tomwalters@0: /*
tomwalters@0:     =============================================================
tomwalters@0:     wdf_tl.c
tomwalters@0:     =============================================================
tomwalters@0: 
tomwalters@0:     Wave digital filter (WDF) implementation of the cochlear 
tomwalters@0:     transmission line (TL) network.
tomwalters@0: 
tomwalters@0:     Author        : Christian Giguere
tomwalters@0:     First written : 01st April, 1991
tomwalters@0:     Last edited   : 18th February, 1994
tomwalters@0: 
tomwalters@0:     References:
tomwalters@0:     (1) C.Giguere and P.C.Woodland (1994). JASA 95(1): 331-342.
tomwalters@0:     (2) A.Fettweis (1986). Proceedings IEEE 74(2): 270-327.
tomwalters@0:     =============================================================
tomwalters@0: */
tomwalters@0: 
tomwalters@0: /***** includes *****/
tomwalters@0: 
tomwalters@0: #include <math.h>
tomwalters@0: #include <stdio.h>
tomwalters@0: #include "stitch.h"
tomwalters@0: #include "source.h"
tomwalters@0: #include "calc.h"
tomwalters@0: #include "calc_tl.h"
tomwalters@0: #include "bank_tl.h" 
tomwalters@0: #include "wdf_tl.h"
tomwalters@0: #include "wdf_ear.h"
tomwalters@0:      
tomwalters@0: /***** defines *****/
tomwalters@0: 
tomwalters@0: #if 0
tomwalters@0: #define   _DEBUG_
tomwalters@0: #define   _OVERFLOW_
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: #if 0
tomwalters@0: #define   _IMPULSE_            /* bypasses input wave with a delta impulse */
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: #if 0			       /* bypasses input wave with a pure tone and dumps  */
tomwalters@0: #define   _BMCURVES_           /* the rms value of basilar membrane velocity on stdout */
tomwalters@0: #endif			       /* (see below for other parameters to set) */
tomwalters@0: 
tomwalters@0: #if 0
tomwalters@0: #define   _EAR_                /* dumps stapes velocity and eardrum pressure on stdout */
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: #define   NUMBER_OF_STATES     (           14 )           /* per TLF segment */
tomwalters@0: #define   NUMBER_OF_COEFFS     (            0 )           /* per TLF segment */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /***** functions *****/
tomwalters@0: 
tomwalters@0: /******************************** WDF set-up function *************************************
tomwalters@0: * name:                 function:
tomwalters@0: *
tomwalters@0: * WDFilterState()       Set-up function for the WDF implementation of the cochlear network 
tomwalters@0: *                       (Giguere and Woodland, 1994, Figs. 3 and 10). It generates the 
tomwalters@0: *                       multiplier coefficients and initializes the states of the filter.
tomwalters@0: *                       It is called by function ``TLF_GenBank()'' in file ``bank_tl.c'' 
tomwalters@0: *                       once for each segment of the transmission line starting from the most 
tomwalters@0: *                       apical segment. It returns a pointer to a structure containing all 
tomwalters@0: *                       the relevant information for the computation of the current filter 
tomwalters@0: *                       segment.
tomwalters@0: *
tomwalters@0: *                       Note: The code that sets up the WDF implementation of the outer ear 
tomwalters@0: *                             and middle ear is located in file ``wdf_ear.c''.
tomwalters@0: ******************************************************************************************/
tomwalters@0: 
tomwalters@0: /************************************ WDFilterState() ************************************/
tomwalters@0: 
tomwalters@0: WDFilterState *WDFilter( samplerate, center_frequency, scale_vel, scale_disp, rho, scala_area, 
tomwalters@0: 			 width, qfactor, mass_per_area, seglength, Lt, warping, active, zov, 
tomwalters@0: 			 OHC_gain, OHC_sat )
tomwalters@0: 
tomwalters@0: double  samplerate, center_frequency, scale_vel, scale_disp ;
tomwalters@0: double  rho, scala_area, width, qfactor, mass_per_area, seglength, Lt ;
tomwalters@0: int     warping, active ;
tomwalters@0: double  *zov, OHC_gain, OHC_sat ;
tomwalters@0: {
tomwalters@0:     static int  apex = 1 ;
tomwalters@0:     DeclareNew( WDFilterState *, filter_state ) ;
tomwalters@0:     CoeffType *cf ;
tomwalters@0:     double  an, bn, mn, d, fs, fn ; 
tomwalters@0:     double  Lsn, Ln, Cn, Rn ;
tomwalters@0:     double  r1, r2, r3, r4, g1, g2, g3, zn, zcn ;
tomwalters@0: 
tomwalters@0:    /***** cochlear parameters for current BM segment *****/
tomwalters@0:     an = scala_area ;                       /* cross-sectional area of each scala (cm2) */
tomwalters@0:     bn = width ;                            /* BM width (cm) */
tomwalters@0:     mn = mass_per_area ;                    /* transversal mass-per-area-of-BM (g/cm2) */
tomwalters@0:     d = seglength ;                         /* BM segment length (cm) */ 
tomwalters@0: 
tomwalters@0:    /***** frequency warping compensation due to bilinear transform *****/
tomwalters@0:     fs = samplerate ;
tomwalters@0:     if ( warping == 0 )
tomwalters@0: 	 fn = center_frequency ;
tomwalters@0:     else
tomwalters@0: 	 fn = fs / Pi * tan( Pi * center_frequency / fs ) ;
tomwalters@0: 
tomwalters@0:    /***** electrical circuit elements (CGS units) for current segment *****/
tomwalters@0:     Ln = mn / ( bn * d ) ;                                      /* shunt inductor */
tomwalters@0:     Cn = 1. / ( ( TwoPi * TwoPi * fn * fn ) * Ln ) ;            /* shunt capacitor */
tomwalters@0:     Rn = sqrt( Ln / Cn ) / qfactor ;                            /* shunt resistor */
tomwalters@0:     Lsn =  ( 2. * rho * d ) / an ;                              /* series inductor */
tomwalters@0: 
tomwalters@0: 
tomwalters@0:    /***** WDF multiplier coefficients for current TL segment *****/
tomwalters@0: 
tomwalters@0:    /*** initialise filter coeffs ***/
tomwalters@0:     filter_state->coeffs = ( char * ) NewZeroedArray( CoeffType, NUMBER_OF_COEFFS, 
tomwalters@0: 						      "wdf_tl.c for filter coeffs\n" ) ;
tomwalters@0:     cf = ( CoeffType * ) filter_state->coeffs ;
tomwalters@0:     
tomwalters@0:     /* Adaptor 25 */
tomwalters@0:     r1 = zcn = 1. / (2. * fs * Cn ) ;                           /* Zcn */
tomwalters@0:     r2 = Rn ;                                                   /* Zrn */
tomwalters@0:     r3 = 2. * fs * Ln ;                                         /* Zln */
tomwalters@0:     r4 = zn = r1 + r2 + r3 ;
tomwalters@0:     filter_state->gamma251 = r1 / r4 ;
tomwalters@0:     filter_state->gamma252 = r2 / r4 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 24 */
tomwalters@0:     g1 = 1. / zn ;
tomwalters@0:     if( apex ) {
tomwalters@0: 	g2 = 1. / ( 2. * fs * Lt ) ;                            /* (1 / Zlt) */
tomwalters@0: 	apex = 0 ; 
tomwalters@0:     }
tomwalters@0:     else
tomwalters@0: 	g2 = 1. / *zov ;                                        /* Zov (at the base) */     
tomwalters@0:     g3 = g1 + g2 ;
tomwalters@0:     filter_state->gamma241 = g1 / g3 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 23 */
tomwalters@0:     r1 = 2. * fs * Lsn ;                                        /* Zlsn */
tomwalters@0:     r2 = 1. / g3 ;
tomwalters@0:     r3 = *zov = r1 + r2 ;                                       /* Zov */
tomwalters@0:     filter_state->gamma231 = r1 / r3 ; 
tomwalters@0: 
tomwalters@0: #ifdef _DEBUG_
tomwalters@0:     fprintf( stderr, "gamma251=%f gamma252=%f\n", filter_state->gamma251
tomwalters@0: 						, filter_state->gamma252 ) ;
tomwalters@0:     fprintf( stderr, "gamma241=%f            \n", filter_state->gamma241 ) ; 
tomwalters@0:     fprintf( stderr, "gamma231=%f            \n", filter_state->gamma231 ) ;
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0:    /***** scaling to convert output to transversal motion of BM segment *****/ 
tomwalters@0:     filter_state->out_scale_disp = scale_disp * Cn / ( 2. * bn * d ) ;
tomwalters@0:     filter_state->out_scale_vel = scale_vel / ( 2. * zcn * bn * d ) ;
tomwalters@0:    
tomwalters@0:    /***** multiplier coefficients for OHC source *****/
tomwalters@0:     filter_state->OHC_sat = OHC_sat * scale_disp / filter_state->out_scale_disp ;
tomwalters@0:     filter_state->OHC_gain = - OHC_gain * Rn * filter_state->OHC_sat / ( 2. * zcn ) ;
tomwalters@0: 
tomwalters@0: #ifdef _BMCURVES_
tomwalters@0:    /***** initialise rms detection *****/
tomwalters@0:     filter_state->rms = 0.0 ;
tomwalters@0:     filter_state->sample = 0 ;
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0:    /***** initialise filter states *****/
tomwalters@0:     filter_state->states = ( char * ) NewZeroedArray( StateType, NUMBER_OF_STATES, 
tomwalters@0: 						      "wdf_tl.c for filter states\n" ) ;
tomwalters@0:     
tomwalters@0:    /***** is channel active for display ? *****/
tomwalters@0:     filter_state->active = active ;
tomwalters@0: 
tomwalters@0:    /***** return *****/
tomwalters@0:     return( filter_state ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*********************************** WDF procedures ***************************************
tomwalters@0: * name:                 function:
tomwalters@0: *
tomwalters@0: * DoWDFdataArray()      WDF realization of the outer ear, middle ear and cochlear networks. 
tomwalters@0: *                       It is called by function ``tlf_bank_callback()'' in file ``bank_tl.c'' 
tomwalters@0: *                       once for a specified number of input points. It computes the BM 
tomwalters@0: *                       displacement or velocity for each segment and keeps track of the 
tomwalters@0: *                       filter states. 
tomwalters@0: *
tomwalters@0: * CloseWDF()            Deletes all private data structures and arrays of the cochlear 
tomwalters@0: *                       transmission line filter upon completion of filtering. It is called 
tomwalters@0: *                       by function ``tlf_bank_callbank()'' in file ``bank_tl.c''.
tomwalters@0: ******************************************************************************************/
tomwalters@0: 
tomwalters@0: /********************************** DoWDFdataArray() **********************************/
tomwalters@0: 
tomwalters@0: void DoWDFdataArray( bankInfo, filter_states, input, output, points, channels,
tomwalters@0: 		     ws, eartube_states, ms, tube_segments )
tomwalters@0: register TLF_BankInfo     *bankInfo ;
tomwalters@0: register WDFilterState    **filter_states ;
tomwalters@0: register DataType         *input ;
tomwalters@0: register DataType         *output ;
tomwalters@0: register int              points, channels ;
tomwalters@0: register WaveWDFstate       *ws ;
tomwalters@0: register EartubeWDFstate    **eartube_states ;
tomwalters@0: register EarmiddleWDFstate  *ms ;
tomwalters@0: register int              tube_segments ;
tomwalters@0: {
tomwalters@0:     register WDFilterState  *fs ;
tomwalters@0:     register EartubeWDFstate *ts ;
tomwalters@0:     register StateType   *st ;   
tomwalters@0:     register StateType   dn, in, out ; 
tomwalters@0:     register StateType   a0, an_1, b1, b2, b3, b63, b161, b202, b212 ;
tomwalters@0:     static   StateType   bn = 0. ;
tomwalters@0:     register CoeffType   *cf ;
tomwalters@0:     register int         decimateCount = bankInfo->decimateCount ;
tomwalters@0:     register int         decimate_factor = bankInfo->decimate_factor ;
tomwalters@0:     register int         channelActive, pointActive, motion = bankInfo->output_var ;
tomwalters@0:     register int         channel = -1, tube_segment = -1 ;
tomwalters@0: #ifdef _BMCURVES_
tomwalters@0:     static   long        sample = -1 ;
tomwalters@0: 	     int         period = 48 ;                        /* in samples */
tomwalters@0: 	     double      samplerate = 48000. ;                /* samples per sec */
tomwalters@0: 	     double      rise_time = 2.0 ;                    /* sec */
tomwalters@0: 	     double      start_time = 4.0 ;                   /* sec */
tomwalters@0: #endif
tomwalters@0: #ifdef _IMPULSE_
tomwalters@0:     static   int         impulse = 10000 ;
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0:     /***** update stapedial muscle state *****/
tomwalters@0: 
tomwalters@0:      ( void ) ms->update_proc( ms ) ; 
tomwalters@0: 
tomwalters@0:     /***** for all points ******/
tomwalters@0: 
tomwalters@0: #ifdef _DEBUG_     
tomwalters@0:      fprintf( stderr, "DoWDFdataArray() = %d points X %d channels\n", points, channels ) ;
tomwalters@0: #endif     
tomwalters@0: 
tomwalters@0:      while( points-- > 0 ) {
tomwalters@0: 
tomwalters@0: 	if( decimateCount == decimate_factor )
tomwalters@0: 	   decimateCount = 0 ;
tomwalters@0: 	pointActive = !decimateCount ;
tomwalters@0: 	decimateCount++ ;
tomwalters@0: 
tomwalters@0: 	output += bankInfo->output_chans * pointActive ; 
tomwalters@0: 
tomwalters@0:        /*** for all TL channels in BACKWARD direction for current input point ***/
tomwalters@0: 
tomwalters@0: 	while( ++channel < channels ) {
tomwalters@0:   
tomwalters@0: 	  /* get states and coefficients */
tomwalters@0: 	  fs = *filter_states++ ;   
tomwalters@0: 	  st = ( StateType * ) fs->states ;
tomwalters@0: 	  cf = ( CoeffType * ) fs->coeffs ; 
tomwalters@0: 
tomwalters@0: 	  /* Adaptor 25 */
tomwalters@0: 	  st[10] = st[1] ;                                    /* a251 */
tomwalters@0: 	  st[3] = st[2] - st[4] - st[1] ;                     /* b254 = a221 */
tomwalters@0: 
tomwalters@0: 	  /* Adaptor 24 */ 
tomwalters@0: 	  st[9] = bn ;                                        /* -a242 */
tomwalters@0: 	  st[5] = fs->gamma241 * ( st[3] + st[9] ) ;          /* b240 */
tomwalters@0: 	  st[6] = st[5] - st[9] ;                             /* b243 = a232 */
tomwalters@0: 
tomwalters@0: 	  /* Adaptor 23 */
tomwalters@0: 	  st[7] = st[8] - st[6] ;                             /* b233 */
tomwalters@0: 	  bn = st[7] ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0:      /*** input  point ***/
tomwalters@0: 
tomwalters@0: 	in = ( StateType ) *input++ ;  
tomwalters@0: 
tomwalters@0: #ifdef _IMPULSE_
tomwalters@0: 	if( impulse != 0 )  {
tomwalters@0: 	   in = ( StateType ) impulse ;
tomwalters@0: 	   impulse = 0 ;
tomwalters@0: 	}
tomwalters@0: 	else 
tomwalters@0: 	   in = ( StateType ) 0 ;
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: #ifdef _BMCURVES_
tomwalters@0: 	if( ++sample < ( long ) ( samplerate * rise_time ) )
tomwalters@0: 	    in = 0.28284 * ( double ) sample / ( samplerate * rise_time ) * 
tomwalters@0: 		 sin( TwoPi * ( double ) ( sample % period ) / ( double ) period ) ;
tomwalters@0: 	else
tomwalters@0: 	    in = 0.28284 * sin( TwoPi * ( double ) ( sample % period ) / ( double ) period ) ;
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0:      /*** computation of outer and middle ear in FORWARD direction for current input point ***/
tomwalters@0: 
tomwalters@0:      /* freefield - external ear */
tomwalters@0: 
tomwalters@0: 	/* get states */
tomwalters@0: 	st = ( StateType * ) ws->states ;      
tomwalters@0: 
tomwalters@0: 	/* Adaptor 0 */
tomwalters@0: 	st[4] = ws->gamma01 * ( st[3] - in ) ;                /* b00  */
tomwalters@0: 	st[5] = st[4] + in + in ;                             /* b03 = -a11 */ 
tomwalters@0: 
tomwalters@0: 	/* Adaptor 2 */
tomwalters@0: 	st[1] = -ws->gamma21 * st[0] ;                        /* b20 */
tomwalters@0: 	st[2] = st[1] ;                                       /* b23 = a12 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 1 */
tomwalters@0: 	st[6] = st[5] - st[2] ;                               /* b13 = a32 */
tomwalters@0: 	an_1 = st[6] ;
tomwalters@0: 
tomwalters@0:      /* concha and auditory canal */
tomwalters@0: 
tomwalters@0: 	while( ++tube_segment < tube_segments ) {
tomwalters@0: 
tomwalters@0: 	   /* get states */
tomwalters@0: 	   ts = *eartube_states++ ;  
tomwalters@0: 	   st = ( StateType * ) ts->states ;
tomwalters@0:  
tomwalters@0: 	   /* Adaptor 3 */
tomwalters@0: 	   st[3] = st[4] - an_1 ;                             /* b33 = a41 */
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 5 */
tomwalters@0: 	   st[0] = -ts->gamma51 * st[1] ;                     /* b50 */
tomwalters@0: 	   st[2] = st[0] + st[1] ;                            /* b53 = a42 */
tomwalters@0:  
tomwalters@0: 	   /* Adaptor 4 */
tomwalters@0: 	   st[5] = st[3] - st[2] ;                            /* -a42+a41 */
tomwalters@0: 	   st[6] = ts->gamma41 * st[5] ;                      /* b40 */
tomwalters@0: 	   st[7] = st[6] + st[2] ;                            /* b43 = a61 */
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 6 */
tomwalters@0: 	   st[8] = st[9] - st[7] ;                            /* b63 = a32 */
tomwalters@0: 	   an_1 = b63 = st[8] ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0:      /* middle ear */
tomwalters@0: 
tomwalters@0: 	/* get states */
tomwalters@0: 	st = ( StateType * ) ms->states ;      
tomwalters@0: 
tomwalters@0: 	/* Adaptor 9 */
tomwalters@0: 	st[31] = st[30] - st[29] ;                           /* b94 = a81 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 8 */
tomwalters@0: 	st[32] = st[35] - st[31] ;                           /* a83-a81 */
tomwalters@0: 	st[33] = -ms->gamma81 * st[32] - ms->gamma82 * st[35] ;   /* b80 */ 
tomwalters@0: 	st[34] = st[33] + st[35] ;                           /* b84 = a71 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 7 */
tomwalters@0: 	st[36] = -st[34] - an_1 ;                            /* b73 = a102 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 13 */
tomwalters@0: 	st[15] = -st[14] ;                                   /* b133 = a122 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 12 */
tomwalters@0: 	st[16] = st[15] + st[18] ;                           /* a122-a121 */
tomwalters@0: 	st[17] = -ms->gamma121 * st[16] ;                    /* b120 */
tomwalters@0: 	st[19] = st[17] + st[15] ;                           /* b123 = a112 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 14 */
tomwalters@0: 	st[21] = -ms->gamma141 * st[20] ;                    /* b140 */
tomwalters@0: 	st[22] = st[21] + st[20] ;                           /* b143 = a111 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 15 */
tomwalters@0: 	st[23] = -st[24] ;                                   /* b153 = a113 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 11 */
tomwalters@0: 	st[25] = -st[22] - st[19] - st[23] ;                 /* b114 = a101 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 10 */
tomwalters@0: 	st[26] = st[25] - st[36] ;                           /* -a102+a101 */
tomwalters@0: 	st[27] = ms->gamma101 * st[26] ;                     /* b100 */
tomwalters@0: 	st[28] = st[27] + st[36] ;                           /* b103 = a161 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 17 */
tomwalters@0: 	st[10] = st[12] - st[11] ;                           /* b174 = a162 */
tomwalters@0: 	
tomwalters@0: 	/* Adaptor 16 */
tomwalters@0: 	st[13] = -st[28] - st[10] ;                          /* b163 = a181 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 19 */
tomwalters@0: 	st[5] = -st[6] ;                                     /* b193 = a182 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 18 */
tomwalters@0: 	st[7] = st[13] - st[5] ;                             /* -a182+a181 */
tomwalters@0: 	st[8] = ms->gamma181 * st[7] ;                       /* b180 */
tomwalters@0: 	st[9] = st[8] + st[5] ;                              /* b183 = a203 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 20 */
tomwalters@0: 	st[40] =  - bn / ms->ratio ;                         /* a202 */
tomwalters@0: 	st[3] = -st[40] - st[9] ;                            /* b204 = a212 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 21 */
tomwalters@0: 	st[0] = st[2] - st[1] - st[3] ;                      /* b214 = a222 */
tomwalters@0: 
tomwalters@0: 
tomwalters@0:      /*** computation of outer and middle ear in BACKWARD direction for current input point ***/
tomwalters@0: 
tomwalters@0:      /* middle ear */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 22 */
tomwalters@0: 	a0 = st[4] + st[0] ;                                  /* a220 */
tomwalters@0: 	st[4] = st[4] - ms->gamma221 * a0 ;                   /* b221 */
tomwalters@0: 	b2 = -st[4] - a0 ;                                    /* b222 = a214 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 21 */
tomwalters@0: 	a0 = b2 - st[0] ;                                     /* a210 */
tomwalters@0: 	st[1] = st[1] - ms->gamma211 * a0 ;                   /* b211 */
tomwalters@0: 	b212 = st[3] - ms->gamma212 * a0 ;                    /* b212 = a204 */
tomwalters@0: 	st[2] = -st[1] - b212 - b2 ;                          /* b213 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 20 */
tomwalters@0: 	a0 = b212 - st[3] ;                                   /* a200 */
tomwalters@0: 	b202 = st[40] - ms->gamma202 * a0 ;                   /* b202 */
tomwalters@0: 	b3 = ms->gamma201 * a0 - b202 - b212 ;                /* b203 = a183 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 18 */
tomwalters@0: 	b2 = st[8] + b3 ;                                     /* b182 = a193 */
tomwalters@0: 	b1 = b2 - st[7] ;                                     /* b181 = a163 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 19 */
tomwalters@0: 	st[6] = st[6] + ms->gamma191 * ( st[5] - b2 ) ;       /* b191 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 16 */
tomwalters@0: 	b161 = st[28] + ms->gamma161 * ( st[13] - b1 ) ;      /* b161 = a103 */
tomwalters@0: 	b2 = -b161 - b1 ;                                     /* b162 = a174 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 17 */
tomwalters@0: 	a0 = b2 - st[10] ;                                    /* a170 */
tomwalters@0: 	st[11] = st[11] - ms->gamma171 * a0 ;                 /* b171 */
tomwalters@0: 	st[12] = ms->gamma172 * a0 - b2 - st[11] ;            /* b173 */
tomwalters@0:  
tomwalters@0: 	/* Adaptor 10 */
tomwalters@0: 	st[37] = st[27] + b161 ;                              /* b102 = a73 */
tomwalters@0: 	st[38] = st[37] - st[26] ;                            /* b101 = a114 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 11 */
tomwalters@0: 	a0 = st[38] - st[25] ;                                /* a110 */
tomwalters@0: 	b1 = st[22] - ms->gamma111 * a0 ;                     /* b111 = a143 */
tomwalters@0: 	b2 = st[19] - ms->gamma112 * a0 ;                     /* b112 = a123 */
tomwalters@0: 	b3 = -b1 -b2 - st[38] ;                               /* b113 = a153 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 15 */
tomwalters@0: 	st[24] = st[24] + ms->gamma151 * ( st[23] - b3 ) ;    /* b151 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 14 */
tomwalters@0: 	st[20] = st[21] + b1 ;                                /* b142 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 12 */
tomwalters@0: 	b2 = st[17] + b2 ;                                    /* b122 = a133 */
tomwalters@0: 	st[18] = b2 + st[16] ;                                /* b121 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 13 */
tomwalters@0: 	st[14] = st[14] + ms->gamma131 * ( st[15] - b2 ) ;    /* b131 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 7 */
tomwalters@0: 	b1 = st[34] + ms->gamma71 * ( st[36] - st[37] ) ;     /* b71 = a84 */
tomwalters@0: 	st[39] = -b1 - st[37] ;                               /* b72 = a63 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 8 */
tomwalters@0: 	st[35] = st[33] + b1 ;                                /* b83 */
tomwalters@0: 	b1 = st[35] + st[32] ;                                /* b81 = a94 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 9 */
tomwalters@0: 	a0 = b1 - st[31] ;                                    /* a90 */
tomwalters@0: 	st[29] = st[29] - ms->gamma91 * a0 ;                  /* b91 */
tomwalters@0: 	st[30] = -st[29] + ms->gamma92 * a0 - b1 ;            /* b93 */
tomwalters@0: 
tomwalters@0:      bn = st[39] ;
tomwalters@0: 
tomwalters@0:      /* concha and auditory canal */
tomwalters@0: 
tomwalters@0: 	while( tube_segment-- > 0 ) {
tomwalters@0: 
tomwalters@0: 	   /* get states */
tomwalters@0: 	   ts = *--eartube_states ;   
tomwalters@0: 	   st = ( StateType * ) ts->states ;
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 6 */
tomwalters@0: 	   b1 = st[7] - ts->gamma61 * ( bn - st[8] ) ;        /* b61 = a43 */
tomwalters@0: 	   st[9] = -b1 - bn ;                                 /* b62 */
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 4 */
tomwalters@0: 	   b2 = st[6] + b1 ;                                  /* b42 = a53 */
tomwalters@0: 	   b1 = b2 - st[5] ;                                  /* b41 = a33 */
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 5 */
tomwalters@0: 	   st[1] = st[0] + b2 ;                               /* b52 */
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 3 */
tomwalters@0: 	   st[4] = ts->gamma31 * ( st[3] - b1 ) - st[4] ;     /* b31 */
tomwalters@0: 	   bn = -st[4] - b1 ;                                 /* b32 = a63 */
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0:      /* freefield - external ear */
tomwalters@0: 
tomwalters@0: 	/* get states */
tomwalters@0: 	st = ( StateType * ) ws->states ;      
tomwalters@0: 
tomwalters@0: 	/* Adaptor 1 */
tomwalters@0: 	b1 = ws->gamma11 * ( st[6] - bn ) - st[5] ;           /* b11 = -a03 */
tomwalters@0: 	b2 = -b1 - bn ;                                       /* b12 = a23 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 2 */
tomwalters@0: 	st[0] = st[1] + b2 + st[0] ;                          /* b61 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 0 */
tomwalters@0: 	st[3] = b1 - st[4] + st[3] ;                          /* -b01 + p */
tomwalters@0: 
tomwalters@0:      an_1 = - b202 * ms->ratio ;    
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 	/* for all channels in FORWARD direction for current input point */
tomwalters@0: 
tomwalters@0: 	 while( channel-- > 0 ) {
tomwalters@0: 
tomwalters@0: 	   /* get states and coefficients */
tomwalters@0: 	   fs = *--filter_states ;   
tomwalters@0: 	   st = ( StateType * ) fs->states ;
tomwalters@0: 	   cf = ( CoeffType * ) fs->coeffs ;
tomwalters@0: 	   channelActive = fs->active ;
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 23 */
tomwalters@0: 	   st[8] = fs->gamma231 * ( st[7] - an_1 ) - st[8] ;    /* b231 */
tomwalters@0: 	   b2 = an_1 + st[8] ;                                  /* -b232 = -a243 */
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 24 */
tomwalters@0: 	   an_1 = b2 - st[5] ;                                  /* -b242 */
tomwalters@0: 	   b1 = -an_1 - st[3] - st[9] ;                         /* b241 = a254 */
tomwalters@0: 	  
tomwalters@0: 	   /* Adaptor 25 */
tomwalters@0: 	   a0 = b1 - st[3] ;                                    /* a250 */
tomwalters@0: 	   st[1] = st[10] - fs->gamma251 * a0 ;                 /* b251 */
tomwalters@0: 	   st[2] = fs->gamma252 * a0 - st[4] - st[1] - b1 ;     /* b253 */
tomwalters@0: 	   dn = st[1] + st[10] ; 
tomwalters@0: 	   in = st[1] - st[10] ;
tomwalters@0: 
tomwalters@0: 	   /* output storage */
tomwalters@0: 	   if( channelActive * pointActive ) {
tomwalters@0: 
tomwalters@0: 	      if( motion == DISPLACEMENT ) 
tomwalters@0: 		 out = fs->out_scale_disp * dn ;
tomwalters@0: 	      else 
tomwalters@0: 		 out = fs->out_scale_vel * in ;
tomwalters@0: #ifdef _OVERFLOW_
tomwalters@0: 	      if( out > _MaxOutput_ || out < _MinOutput_ )
tomwalters@0: 		 fprintf( stderr, "Overflow error in BMM output=%e\n", ( double ) out ) ;     
tomwalters@0: #endif
tomwalters@0: 	      *(--output) = ( DataType ) ( out ) ; 
tomwalters@0: 	   }
tomwalters@0: 
tomwalters@0: 	   /* OHC nonlinear voltage source */
tomwalters@0: 	      if( dn < 0. )
tomwalters@0: 		 dn = -dn ;
tomwalters@0: 	      st[4] = fs->OHC_gain * in / ( fs->OHC_sat + dn ) ;  /* -Vn(ohc) */
tomwalters@0: 
tomwalters@0: #ifdef _BMCURVES_
tomwalters@0: 	   if( sample > ( long ) ( start_time * samplerate ) ) { 
tomwalters@0: 	      in = fs->out_scale_vel * in ;                     
tomwalters@0: 	      fs->rms = fs->rms + in*in ; 
tomwalters@0: 	      fs->sample++ ;
tomwalters@0: 	   }
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: 	}
tomwalters@0: 	bn = -an_1 ;                                     /* apical boundary condition */
tomwalters@0: 	output += bankInfo->output_chans * pointActive ;
tomwalters@0: 
tomwalters@0: #ifdef _EAR_
tomwalters@0: 	st = ( StateType * ) ms->states ;
tomwalters@0: 	printf( "stapes=%e  eardrum=%e\n", 0.5 * ( st[40] - b202 ), 0.5 * ( st[39] + b63 ) ) ; 
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     return ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /******************************* CloseWDF() *************************************/
tomwalters@0: 
tomwalters@0: void CloseWDF( states, channels, bank )
tomwalters@0:   register WDFilterState      **states ;
tomwalters@0:   register int                channels ;
tomwalters@0:   register TLF_BankInfo       *bank ;
tomwalters@0: {
tomwalters@0:   int channel ;
tomwalters@0: 
tomwalters@0:   for( channel = 0 ; channel < channels ; channel++ ) {
tomwalters@0: #ifdef _BMCURVES_
tomwalters@0:      printf( "%d  %e\n", channel+1, sqrt(states[channel]->rms / states[channel]->sample) ) ;
tomwalters@0: #endif
tomwalters@0:      Delete( states[channel]->states ) ;
tomwalters@0:      Delete( states[channel]->coeffs ) ;
tomwalters@0:   }
tomwalters@0:   Delete( states ) ;
tomwalters@0:   
tomwalters@0:   Delete( bank ) ;
tomwalters@0: 
tomwalters@0:   return ;
tomwalters@0: 
tomwalters@0: }
tomwalters@0: