Mercurial > hg > aim92
diff wdf/wdf_tl.c @ 0:5242703e91d3 tip
Initial checkin for AIM92 aimR8.2 (last updated May 1997).
| author | tomwalters |
|---|---|
| date | Fri, 20 May 2011 15:19:45 +0100 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/wdf/wdf_tl.c Fri May 20 15:19:45 2011 +0100 @@ -0,0 +1,637 @@ +/* + 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 ; + +} +