Mercurial > hg > aim92
comparison 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 |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:5242703e91d3 |
---|---|
1 /* | |
2 Copyright (c) Applied Psychology Unit, Medical Research Council. 1988, 1989 | |
3 =========================================================================== | |
4 | |
5 Permission to use, copy, modify, and distribute this software without fee | |
6 is hereby granted for research purposes, provided that this copyright | |
7 notice appears in all copies and in all supporting documentation, and that | |
8 the software is not redistributed for any fee (except for a nominal shipping | |
9 charge). Anyone wanting to incorporate all or part of this software in a | |
10 commercial product must obtain a license from the Medical Research Council. | |
11 | |
12 The MRC makes no representations about the suitability of this | |
13 software for any purpose. It is provided "as is" without express or implied | |
14 warranty. | |
15 | |
16 THE MRC DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING | |
17 ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL THE | |
18 A.P.U. BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY | |
19 DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN | |
20 AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF | |
21 OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. | |
22 */ | |
23 | |
24 /* | |
25 Acknowledgment: | |
26 ============== | |
27 | |
28 The source code provided in this file was originally developed by | |
29 Christian Giguere as part of a Ph.D degree at the Department of | |
30 Engineering of the University of Cambridge from April 1990 to | |
31 November 1993. The code was subsequently adapted under a grant | |
32 from the Hearing Research Trust for full compatibility with | |
33 AIM Release 6.15. | |
34 | |
35 Christian Giguere 25/03/94 | |
36 | |
37 */ | |
38 | |
39 /* | |
40 ============================================================= | |
41 wdf_tl.c | |
42 ============================================================= | |
43 | |
44 Wave digital filter (WDF) implementation of the cochlear | |
45 transmission line (TL) network. | |
46 | |
47 Author : Christian Giguere | |
48 First written : 01st April, 1991 | |
49 Last edited : 18th February, 1994 | |
50 | |
51 References: | |
52 (1) C.Giguere and P.C.Woodland (1994). JASA 95(1): 331-342. | |
53 (2) A.Fettweis (1986). Proceedings IEEE 74(2): 270-327. | |
54 ============================================================= | |
55 */ | |
56 | |
57 /***** includes *****/ | |
58 | |
59 #include <math.h> | |
60 #include <stdio.h> | |
61 #include "stitch.h" | |
62 #include "source.h" | |
63 #include "calc.h" | |
64 #include "calc_tl.h" | |
65 #include "bank_tl.h" | |
66 #include "wdf_tl.h" | |
67 #include "wdf_ear.h" | |
68 | |
69 /***** defines *****/ | |
70 | |
71 #if 0 | |
72 #define _DEBUG_ | |
73 #define _OVERFLOW_ | |
74 #endif | |
75 | |
76 #if 0 | |
77 #define _IMPULSE_ /* bypasses input wave with a delta impulse */ | |
78 #endif | |
79 | |
80 #if 0 /* bypasses input wave with a pure tone and dumps */ | |
81 #define _BMCURVES_ /* the rms value of basilar membrane velocity on stdout */ | |
82 #endif /* (see below for other parameters to set) */ | |
83 | |
84 #if 0 | |
85 #define _EAR_ /* dumps stapes velocity and eardrum pressure on stdout */ | |
86 #endif | |
87 | |
88 #define NUMBER_OF_STATES ( 14 ) /* per TLF segment */ | |
89 #define NUMBER_OF_COEFFS ( 0 ) /* per TLF segment */ | |
90 | |
91 | |
92 /***** functions *****/ | |
93 | |
94 /******************************** WDF set-up function ************************************* | |
95 * name: function: | |
96 * | |
97 * WDFilterState() Set-up function for the WDF implementation of the cochlear network | |
98 * (Giguere and Woodland, 1994, Figs. 3 and 10). It generates the | |
99 * multiplier coefficients and initializes the states of the filter. | |
100 * It is called by function ``TLF_GenBank()'' in file ``bank_tl.c'' | |
101 * once for each segment of the transmission line starting from the most | |
102 * apical segment. It returns a pointer to a structure containing all | |
103 * the relevant information for the computation of the current filter | |
104 * segment. | |
105 * | |
106 * Note: The code that sets up the WDF implementation of the outer ear | |
107 * and middle ear is located in file ``wdf_ear.c''. | |
108 ******************************************************************************************/ | |
109 | |
110 /************************************ WDFilterState() ************************************/ | |
111 | |
112 WDFilterState *WDFilter( samplerate, center_frequency, scale_vel, scale_disp, rho, scala_area, | |
113 width, qfactor, mass_per_area, seglength, Lt, warping, active, zov, | |
114 OHC_gain, OHC_sat ) | |
115 | |
116 double samplerate, center_frequency, scale_vel, scale_disp ; | |
117 double rho, scala_area, width, qfactor, mass_per_area, seglength, Lt ; | |
118 int warping, active ; | |
119 double *zov, OHC_gain, OHC_sat ; | |
120 { | |
121 static int apex = 1 ; | |
122 DeclareNew( WDFilterState *, filter_state ) ; | |
123 CoeffType *cf ; | |
124 double an, bn, mn, d, fs, fn ; | |
125 double Lsn, Ln, Cn, Rn ; | |
126 double r1, r2, r3, r4, g1, g2, g3, zn, zcn ; | |
127 | |
128 /***** cochlear parameters for current BM segment *****/ | |
129 an = scala_area ; /* cross-sectional area of each scala (cm2) */ | |
130 bn = width ; /* BM width (cm) */ | |
131 mn = mass_per_area ; /* transversal mass-per-area-of-BM (g/cm2) */ | |
132 d = seglength ; /* BM segment length (cm) */ | |
133 | |
134 /***** frequency warping compensation due to bilinear transform *****/ | |
135 fs = samplerate ; | |
136 if ( warping == 0 ) | |
137 fn = center_frequency ; | |
138 else | |
139 fn = fs / Pi * tan( Pi * center_frequency / fs ) ; | |
140 | |
141 /***** electrical circuit elements (CGS units) for current segment *****/ | |
142 Ln = mn / ( bn * d ) ; /* shunt inductor */ | |
143 Cn = 1. / ( ( TwoPi * TwoPi * fn * fn ) * Ln ) ; /* shunt capacitor */ | |
144 Rn = sqrt( Ln / Cn ) / qfactor ; /* shunt resistor */ | |
145 Lsn = ( 2. * rho * d ) / an ; /* series inductor */ | |
146 | |
147 | |
148 /***** WDF multiplier coefficients for current TL segment *****/ | |
149 | |
150 /*** initialise filter coeffs ***/ | |
151 filter_state->coeffs = ( char * ) NewZeroedArray( CoeffType, NUMBER_OF_COEFFS, | |
152 "wdf_tl.c for filter coeffs\n" ) ; | |
153 cf = ( CoeffType * ) filter_state->coeffs ; | |
154 | |
155 /* Adaptor 25 */ | |
156 r1 = zcn = 1. / (2. * fs * Cn ) ; /* Zcn */ | |
157 r2 = Rn ; /* Zrn */ | |
158 r3 = 2. * fs * Ln ; /* Zln */ | |
159 r4 = zn = r1 + r2 + r3 ; | |
160 filter_state->gamma251 = r1 / r4 ; | |
161 filter_state->gamma252 = r2 / r4 ; | |
162 | |
163 /* Adaptor 24 */ | |
164 g1 = 1. / zn ; | |
165 if( apex ) { | |
166 g2 = 1. / ( 2. * fs * Lt ) ; /* (1 / Zlt) */ | |
167 apex = 0 ; | |
168 } | |
169 else | |
170 g2 = 1. / *zov ; /* Zov (at the base) */ | |
171 g3 = g1 + g2 ; | |
172 filter_state->gamma241 = g1 / g3 ; | |
173 | |
174 /* Adaptor 23 */ | |
175 r1 = 2. * fs * Lsn ; /* Zlsn */ | |
176 r2 = 1. / g3 ; | |
177 r3 = *zov = r1 + r2 ; /* Zov */ | |
178 filter_state->gamma231 = r1 / r3 ; | |
179 | |
180 #ifdef _DEBUG_ | |
181 fprintf( stderr, "gamma251=%f gamma252=%f\n", filter_state->gamma251 | |
182 , filter_state->gamma252 ) ; | |
183 fprintf( stderr, "gamma241=%f \n", filter_state->gamma241 ) ; | |
184 fprintf( stderr, "gamma231=%f \n", filter_state->gamma231 ) ; | |
185 #endif | |
186 | |
187 /***** scaling to convert output to transversal motion of BM segment *****/ | |
188 filter_state->out_scale_disp = scale_disp * Cn / ( 2. * bn * d ) ; | |
189 filter_state->out_scale_vel = scale_vel / ( 2. * zcn * bn * d ) ; | |
190 | |
191 /***** multiplier coefficients for OHC source *****/ | |
192 filter_state->OHC_sat = OHC_sat * scale_disp / filter_state->out_scale_disp ; | |
193 filter_state->OHC_gain = - OHC_gain * Rn * filter_state->OHC_sat / ( 2. * zcn ) ; | |
194 | |
195 #ifdef _BMCURVES_ | |
196 /***** initialise rms detection *****/ | |
197 filter_state->rms = 0.0 ; | |
198 filter_state->sample = 0 ; | |
199 #endif | |
200 | |
201 /***** initialise filter states *****/ | |
202 filter_state->states = ( char * ) NewZeroedArray( StateType, NUMBER_OF_STATES, | |
203 "wdf_tl.c for filter states\n" ) ; | |
204 | |
205 /***** is channel active for display ? *****/ | |
206 filter_state->active = active ; | |
207 | |
208 /***** return *****/ | |
209 return( filter_state ) ; | |
210 } | |
211 | |
212 | |
213 /*********************************** WDF procedures *************************************** | |
214 * name: function: | |
215 * | |
216 * DoWDFdataArray() WDF realization of the outer ear, middle ear and cochlear networks. | |
217 * It is called by function ``tlf_bank_callback()'' in file ``bank_tl.c'' | |
218 * once for a specified number of input points. It computes the BM | |
219 * displacement or velocity for each segment and keeps track of the | |
220 * filter states. | |
221 * | |
222 * CloseWDF() Deletes all private data structures and arrays of the cochlear | |
223 * transmission line filter upon completion of filtering. It is called | |
224 * by function ``tlf_bank_callbank()'' in file ``bank_tl.c''. | |
225 ******************************************************************************************/ | |
226 | |
227 /********************************** DoWDFdataArray() **********************************/ | |
228 | |
229 void DoWDFdataArray( bankInfo, filter_states, input, output, points, channels, | |
230 ws, eartube_states, ms, tube_segments ) | |
231 register TLF_BankInfo *bankInfo ; | |
232 register WDFilterState **filter_states ; | |
233 register DataType *input ; | |
234 register DataType *output ; | |
235 register int points, channels ; | |
236 register WaveWDFstate *ws ; | |
237 register EartubeWDFstate **eartube_states ; | |
238 register EarmiddleWDFstate *ms ; | |
239 register int tube_segments ; | |
240 { | |
241 register WDFilterState *fs ; | |
242 register EartubeWDFstate *ts ; | |
243 register StateType *st ; | |
244 register StateType dn, in, out ; | |
245 register StateType a0, an_1, b1, b2, b3, b63, b161, b202, b212 ; | |
246 static StateType bn = 0. ; | |
247 register CoeffType *cf ; | |
248 register int decimateCount = bankInfo->decimateCount ; | |
249 register int decimate_factor = bankInfo->decimate_factor ; | |
250 register int channelActive, pointActive, motion = bankInfo->output_var ; | |
251 register int channel = -1, tube_segment = -1 ; | |
252 #ifdef _BMCURVES_ | |
253 static long sample = -1 ; | |
254 int period = 48 ; /* in samples */ | |
255 double samplerate = 48000. ; /* samples per sec */ | |
256 double rise_time = 2.0 ; /* sec */ | |
257 double start_time = 4.0 ; /* sec */ | |
258 #endif | |
259 #ifdef _IMPULSE_ | |
260 static int impulse = 10000 ; | |
261 #endif | |
262 | |
263 /***** update stapedial muscle state *****/ | |
264 | |
265 ( void ) ms->update_proc( ms ) ; | |
266 | |
267 /***** for all points ******/ | |
268 | |
269 #ifdef _DEBUG_ | |
270 fprintf( stderr, "DoWDFdataArray() = %d points X %d channels\n", points, channels ) ; | |
271 #endif | |
272 | |
273 while( points-- > 0 ) { | |
274 | |
275 if( decimateCount == decimate_factor ) | |
276 decimateCount = 0 ; | |
277 pointActive = !decimateCount ; | |
278 decimateCount++ ; | |
279 | |
280 output += bankInfo->output_chans * pointActive ; | |
281 | |
282 /*** for all TL channels in BACKWARD direction for current input point ***/ | |
283 | |
284 while( ++channel < channels ) { | |
285 | |
286 /* get states and coefficients */ | |
287 fs = *filter_states++ ; | |
288 st = ( StateType * ) fs->states ; | |
289 cf = ( CoeffType * ) fs->coeffs ; | |
290 | |
291 /* Adaptor 25 */ | |
292 st[10] = st[1] ; /* a251 */ | |
293 st[3] = st[2] - st[4] - st[1] ; /* b254 = a221 */ | |
294 | |
295 /* Adaptor 24 */ | |
296 st[9] = bn ; /* -a242 */ | |
297 st[5] = fs->gamma241 * ( st[3] + st[9] ) ; /* b240 */ | |
298 st[6] = st[5] - st[9] ; /* b243 = a232 */ | |
299 | |
300 /* Adaptor 23 */ | |
301 st[7] = st[8] - st[6] ; /* b233 */ | |
302 bn = st[7] ; | |
303 } | |
304 | |
305 /*** input point ***/ | |
306 | |
307 in = ( StateType ) *input++ ; | |
308 | |
309 #ifdef _IMPULSE_ | |
310 if( impulse != 0 ) { | |
311 in = ( StateType ) impulse ; | |
312 impulse = 0 ; | |
313 } | |
314 else | |
315 in = ( StateType ) 0 ; | |
316 #endif | |
317 | |
318 #ifdef _BMCURVES_ | |
319 if( ++sample < ( long ) ( samplerate * rise_time ) ) | |
320 in = 0.28284 * ( double ) sample / ( samplerate * rise_time ) * | |
321 sin( TwoPi * ( double ) ( sample % period ) / ( double ) period ) ; | |
322 else | |
323 in = 0.28284 * sin( TwoPi * ( double ) ( sample % period ) / ( double ) period ) ; | |
324 #endif | |
325 | |
326 /*** computation of outer and middle ear in FORWARD direction for current input point ***/ | |
327 | |
328 /* freefield - external ear */ | |
329 | |
330 /* get states */ | |
331 st = ( StateType * ) ws->states ; | |
332 | |
333 /* Adaptor 0 */ | |
334 st[4] = ws->gamma01 * ( st[3] - in ) ; /* b00 */ | |
335 st[5] = st[4] + in + in ; /* b03 = -a11 */ | |
336 | |
337 /* Adaptor 2 */ | |
338 st[1] = -ws->gamma21 * st[0] ; /* b20 */ | |
339 st[2] = st[1] ; /* b23 = a12 */ | |
340 | |
341 /* Adaptor 1 */ | |
342 st[6] = st[5] - st[2] ; /* b13 = a32 */ | |
343 an_1 = st[6] ; | |
344 | |
345 /* concha and auditory canal */ | |
346 | |
347 while( ++tube_segment < tube_segments ) { | |
348 | |
349 /* get states */ | |
350 ts = *eartube_states++ ; | |
351 st = ( StateType * ) ts->states ; | |
352 | |
353 /* Adaptor 3 */ | |
354 st[3] = st[4] - an_1 ; /* b33 = a41 */ | |
355 | |
356 /* Adaptor 5 */ | |
357 st[0] = -ts->gamma51 * st[1] ; /* b50 */ | |
358 st[2] = st[0] + st[1] ; /* b53 = a42 */ | |
359 | |
360 /* Adaptor 4 */ | |
361 st[5] = st[3] - st[2] ; /* -a42+a41 */ | |
362 st[6] = ts->gamma41 * st[5] ; /* b40 */ | |
363 st[7] = st[6] + st[2] ; /* b43 = a61 */ | |
364 | |
365 /* Adaptor 6 */ | |
366 st[8] = st[9] - st[7] ; /* b63 = a32 */ | |
367 an_1 = b63 = st[8] ; | |
368 } | |
369 | |
370 /* middle ear */ | |
371 | |
372 /* get states */ | |
373 st = ( StateType * ) ms->states ; | |
374 | |
375 /* Adaptor 9 */ | |
376 st[31] = st[30] - st[29] ; /* b94 = a81 */ | |
377 | |
378 /* Adaptor 8 */ | |
379 st[32] = st[35] - st[31] ; /* a83-a81 */ | |
380 st[33] = -ms->gamma81 * st[32] - ms->gamma82 * st[35] ; /* b80 */ | |
381 st[34] = st[33] + st[35] ; /* b84 = a71 */ | |
382 | |
383 /* Adaptor 7 */ | |
384 st[36] = -st[34] - an_1 ; /* b73 = a102 */ | |
385 | |
386 /* Adaptor 13 */ | |
387 st[15] = -st[14] ; /* b133 = a122 */ | |
388 | |
389 /* Adaptor 12 */ | |
390 st[16] = st[15] + st[18] ; /* a122-a121 */ | |
391 st[17] = -ms->gamma121 * st[16] ; /* b120 */ | |
392 st[19] = st[17] + st[15] ; /* b123 = a112 */ | |
393 | |
394 /* Adaptor 14 */ | |
395 st[21] = -ms->gamma141 * st[20] ; /* b140 */ | |
396 st[22] = st[21] + st[20] ; /* b143 = a111 */ | |
397 | |
398 /* Adaptor 15 */ | |
399 st[23] = -st[24] ; /* b153 = a113 */ | |
400 | |
401 /* Adaptor 11 */ | |
402 st[25] = -st[22] - st[19] - st[23] ; /* b114 = a101 */ | |
403 | |
404 /* Adaptor 10 */ | |
405 st[26] = st[25] - st[36] ; /* -a102+a101 */ | |
406 st[27] = ms->gamma101 * st[26] ; /* b100 */ | |
407 st[28] = st[27] + st[36] ; /* b103 = a161 */ | |
408 | |
409 /* Adaptor 17 */ | |
410 st[10] = st[12] - st[11] ; /* b174 = a162 */ | |
411 | |
412 /* Adaptor 16 */ | |
413 st[13] = -st[28] - st[10] ; /* b163 = a181 */ | |
414 | |
415 /* Adaptor 19 */ | |
416 st[5] = -st[6] ; /* b193 = a182 */ | |
417 | |
418 /* Adaptor 18 */ | |
419 st[7] = st[13] - st[5] ; /* -a182+a181 */ | |
420 st[8] = ms->gamma181 * st[7] ; /* b180 */ | |
421 st[9] = st[8] + st[5] ; /* b183 = a203 */ | |
422 | |
423 /* Adaptor 20 */ | |
424 st[40] = - bn / ms->ratio ; /* a202 */ | |
425 st[3] = -st[40] - st[9] ; /* b204 = a212 */ | |
426 | |
427 /* Adaptor 21 */ | |
428 st[0] = st[2] - st[1] - st[3] ; /* b214 = a222 */ | |
429 | |
430 | |
431 /*** computation of outer and middle ear in BACKWARD direction for current input point ***/ | |
432 | |
433 /* middle ear */ | |
434 | |
435 /* Adaptor 22 */ | |
436 a0 = st[4] + st[0] ; /* a220 */ | |
437 st[4] = st[4] - ms->gamma221 * a0 ; /* b221 */ | |
438 b2 = -st[4] - a0 ; /* b222 = a214 */ | |
439 | |
440 /* Adaptor 21 */ | |
441 a0 = b2 - st[0] ; /* a210 */ | |
442 st[1] = st[1] - ms->gamma211 * a0 ; /* b211 */ | |
443 b212 = st[3] - ms->gamma212 * a0 ; /* b212 = a204 */ | |
444 st[2] = -st[1] - b212 - b2 ; /* b213 */ | |
445 | |
446 /* Adaptor 20 */ | |
447 a0 = b212 - st[3] ; /* a200 */ | |
448 b202 = st[40] - ms->gamma202 * a0 ; /* b202 */ | |
449 b3 = ms->gamma201 * a0 - b202 - b212 ; /* b203 = a183 */ | |
450 | |
451 /* Adaptor 18 */ | |
452 b2 = st[8] + b3 ; /* b182 = a193 */ | |
453 b1 = b2 - st[7] ; /* b181 = a163 */ | |
454 | |
455 /* Adaptor 19 */ | |
456 st[6] = st[6] + ms->gamma191 * ( st[5] - b2 ) ; /* b191 */ | |
457 | |
458 /* Adaptor 16 */ | |
459 b161 = st[28] + ms->gamma161 * ( st[13] - b1 ) ; /* b161 = a103 */ | |
460 b2 = -b161 - b1 ; /* b162 = a174 */ | |
461 | |
462 /* Adaptor 17 */ | |
463 a0 = b2 - st[10] ; /* a170 */ | |
464 st[11] = st[11] - ms->gamma171 * a0 ; /* b171 */ | |
465 st[12] = ms->gamma172 * a0 - b2 - st[11] ; /* b173 */ | |
466 | |
467 /* Adaptor 10 */ | |
468 st[37] = st[27] + b161 ; /* b102 = a73 */ | |
469 st[38] = st[37] - st[26] ; /* b101 = a114 */ | |
470 | |
471 /* Adaptor 11 */ | |
472 a0 = st[38] - st[25] ; /* a110 */ | |
473 b1 = st[22] - ms->gamma111 * a0 ; /* b111 = a143 */ | |
474 b2 = st[19] - ms->gamma112 * a0 ; /* b112 = a123 */ | |
475 b3 = -b1 -b2 - st[38] ; /* b113 = a153 */ | |
476 | |
477 /* Adaptor 15 */ | |
478 st[24] = st[24] + ms->gamma151 * ( st[23] - b3 ) ; /* b151 */ | |
479 | |
480 /* Adaptor 14 */ | |
481 st[20] = st[21] + b1 ; /* b142 */ | |
482 | |
483 /* Adaptor 12 */ | |
484 b2 = st[17] + b2 ; /* b122 = a133 */ | |
485 st[18] = b2 + st[16] ; /* b121 */ | |
486 | |
487 /* Adaptor 13 */ | |
488 st[14] = st[14] + ms->gamma131 * ( st[15] - b2 ) ; /* b131 */ | |
489 | |
490 /* Adaptor 7 */ | |
491 b1 = st[34] + ms->gamma71 * ( st[36] - st[37] ) ; /* b71 = a84 */ | |
492 st[39] = -b1 - st[37] ; /* b72 = a63 */ | |
493 | |
494 /* Adaptor 8 */ | |
495 st[35] = st[33] + b1 ; /* b83 */ | |
496 b1 = st[35] + st[32] ; /* b81 = a94 */ | |
497 | |
498 /* Adaptor 9 */ | |
499 a0 = b1 - st[31] ; /* a90 */ | |
500 st[29] = st[29] - ms->gamma91 * a0 ; /* b91 */ | |
501 st[30] = -st[29] + ms->gamma92 * a0 - b1 ; /* b93 */ | |
502 | |
503 bn = st[39] ; | |
504 | |
505 /* concha and auditory canal */ | |
506 | |
507 while( tube_segment-- > 0 ) { | |
508 | |
509 /* get states */ | |
510 ts = *--eartube_states ; | |
511 st = ( StateType * ) ts->states ; | |
512 | |
513 /* Adaptor 6 */ | |
514 b1 = st[7] - ts->gamma61 * ( bn - st[8] ) ; /* b61 = a43 */ | |
515 st[9] = -b1 - bn ; /* b62 */ | |
516 | |
517 /* Adaptor 4 */ | |
518 b2 = st[6] + b1 ; /* b42 = a53 */ | |
519 b1 = b2 - st[5] ; /* b41 = a33 */ | |
520 | |
521 /* Adaptor 5 */ | |
522 st[1] = st[0] + b2 ; /* b52 */ | |
523 | |
524 /* Adaptor 3 */ | |
525 st[4] = ts->gamma31 * ( st[3] - b1 ) - st[4] ; /* b31 */ | |
526 bn = -st[4] - b1 ; /* b32 = a63 */ | |
527 } | |
528 | |
529 /* freefield - external ear */ | |
530 | |
531 /* get states */ | |
532 st = ( StateType * ) ws->states ; | |
533 | |
534 /* Adaptor 1 */ | |
535 b1 = ws->gamma11 * ( st[6] - bn ) - st[5] ; /* b11 = -a03 */ | |
536 b2 = -b1 - bn ; /* b12 = a23 */ | |
537 | |
538 /* Adaptor 2 */ | |
539 st[0] = st[1] + b2 + st[0] ; /* b61 */ | |
540 | |
541 /* Adaptor 0 */ | |
542 st[3] = b1 - st[4] + st[3] ; /* -b01 + p */ | |
543 | |
544 an_1 = - b202 * ms->ratio ; | |
545 | |
546 | |
547 /* for all channels in FORWARD direction for current input point */ | |
548 | |
549 while( channel-- > 0 ) { | |
550 | |
551 /* get states and coefficients */ | |
552 fs = *--filter_states ; | |
553 st = ( StateType * ) fs->states ; | |
554 cf = ( CoeffType * ) fs->coeffs ; | |
555 channelActive = fs->active ; | |
556 | |
557 /* Adaptor 23 */ | |
558 st[8] = fs->gamma231 * ( st[7] - an_1 ) - st[8] ; /* b231 */ | |
559 b2 = an_1 + st[8] ; /* -b232 = -a243 */ | |
560 | |
561 /* Adaptor 24 */ | |
562 an_1 = b2 - st[5] ; /* -b242 */ | |
563 b1 = -an_1 - st[3] - st[9] ; /* b241 = a254 */ | |
564 | |
565 /* Adaptor 25 */ | |
566 a0 = b1 - st[3] ; /* a250 */ | |
567 st[1] = st[10] - fs->gamma251 * a0 ; /* b251 */ | |
568 st[2] = fs->gamma252 * a0 - st[4] - st[1] - b1 ; /* b253 */ | |
569 dn = st[1] + st[10] ; | |
570 in = st[1] - st[10] ; | |
571 | |
572 /* output storage */ | |
573 if( channelActive * pointActive ) { | |
574 | |
575 if( motion == DISPLACEMENT ) | |
576 out = fs->out_scale_disp * dn ; | |
577 else | |
578 out = fs->out_scale_vel * in ; | |
579 #ifdef _OVERFLOW_ | |
580 if( out > _MaxOutput_ || out < _MinOutput_ ) | |
581 fprintf( stderr, "Overflow error in BMM output=%e\n", ( double ) out ) ; | |
582 #endif | |
583 *(--output) = ( DataType ) ( out ) ; | |
584 } | |
585 | |
586 /* OHC nonlinear voltage source */ | |
587 if( dn < 0. ) | |
588 dn = -dn ; | |
589 st[4] = fs->OHC_gain * in / ( fs->OHC_sat + dn ) ; /* -Vn(ohc) */ | |
590 | |
591 #ifdef _BMCURVES_ | |
592 if( sample > ( long ) ( start_time * samplerate ) ) { | |
593 in = fs->out_scale_vel * in ; | |
594 fs->rms = fs->rms + in*in ; | |
595 fs->sample++ ; | |
596 } | |
597 #endif | |
598 | |
599 } | |
600 bn = -an_1 ; /* apical boundary condition */ | |
601 output += bankInfo->output_chans * pointActive ; | |
602 | |
603 #ifdef _EAR_ | |
604 st = ( StateType * ) ms->states ; | |
605 printf( "stapes=%e eardrum=%e\n", 0.5 * ( st[40] - b202 ), 0.5 * ( st[39] + b63 ) ) ; | |
606 #endif | |
607 | |
608 } | |
609 | |
610 return ; | |
611 } | |
612 | |
613 | |
614 /******************************* CloseWDF() *************************************/ | |
615 | |
616 void CloseWDF( states, channels, bank ) | |
617 register WDFilterState **states ; | |
618 register int channels ; | |
619 register TLF_BankInfo *bank ; | |
620 { | |
621 int channel ; | |
622 | |
623 for( channel = 0 ; channel < channels ; channel++ ) { | |
624 #ifdef _BMCURVES_ | |
625 printf( "%d %e\n", channel+1, sqrt(states[channel]->rms / states[channel]->sample) ) ; | |
626 #endif | |
627 Delete( states[channel]->states ) ; | |
628 Delete( states[channel]->coeffs ) ; | |
629 } | |
630 Delete( states ) ; | |
631 | |
632 Delete( bank ) ; | |
633 | |
634 return ; | |
635 | |
636 } | |
637 |