Mercurial > hg > aim92
diff wdf/meddis.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/meddis.c Fri May 20 15:19:45 2011 +0100 @@ -0,0 +1,491 @@ +/* + 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 + +*/ + +/* + =========================================================== + meddis.c + =========================================================== + + Implementation of Ray Meddis haircell model. + + Author : Christian Giguere + First written : 10th September, 1991 + Last edited : 07th March, 1994 + + References: + (1) C.Giguere and P.C.Woodland (1994). JASA 95(1): 331-342. + (2) R.Meddis et al. (1990). JASA 87(4): 1813-1816. + (3) R.Meddis (1988). JASA 83(3): 1056-1063. + + Note - This is a complete re-write of the original file + "haircell.c" Release 3 (20th September 1988) from + J. Holdsworth and the Applied Psychology Unit. + ============================================================ +*/ + + +/***** includes *****/ + +#include <math.h> +#include <stdio.h> +#include "stitch.h" +#include "source.h" +#include "calc.h" +#include "calc_tl.h" +#include "meddis.h" + + +/***** defines *****/ + +#if 0 +#define _DEBUG_ +#define _OVERFLOW_ +#endif + +#if 1 /* define numerical implementation */ +#define _WDF_ /* Giguere and Woodland (1994) */ +#else +#define _FORWARD_DIFFERENCE_ /* Meddis et al. (1990) */ +#endif + +#if 0 +#define _CLAMPING_ /* clamping of reservoirs */ +#endif + +#define NUMBER_OF_COEFFS ( 12 ) +#define NUMBER_OF_STATES ( 12 ) + + +/* Feedback parameter */ + +#define InputGain_max ( 24.0 ) /* max coupling gain in dB */ + + +/* Meddis et al. (1990) model parameters */ + +#define A_medium ( 10.0 ) /* medium-spontaneous rate fibre */ +#define B_medium ( 3000.0 ) +#define g_medium ( 1000.0 ) + +#define A_high ( 5.0 ) /* high-spontaneous rate fibre */ +#define B_high ( 300.0 ) +#define g_high ( 2000.0 ) + +#define y_value ( 5.05 ) /* replenishment rate */ +#define l_value ( 2500 ) /* rate of loss from the cleft */ +#define r_value ( 6580 ) /* rate of return from the cleft */ +#define x_value ( 66.31 ) /* rate of release from w store */ +#define m_value ( 1.0 ) /* max number of quanta */ +#define h_value ( 50000. ) /* firing-rate factor */ + + +/***** private data structures *****/ + +typedef struct _haircell_info HaircellInfo ; +typedef struct _haircell_bank HaircellBank ; + +struct _haircell_info { + int number_of_states ; + double output_gain ; + StateType *states ; + } ; + +struct _haircell_bank { + int channels ; + double *inputGain_ratio ; + double (*update_proc)() ; /* procedure to update cells input gain */ + int number_of_coeffs ; + CoeffType *coeffs ; + HaircellInfo **cells ; + } ; + + +/***** external variables *****/ + +double *InputGain_ratio ; + + +/***** functions *****/ + +/******************************************************************************* +* name: function: +* +* NewHaircells() Set-up function for the implementation of a bank of inner +* haircells based on the model of Ray Meddis. It generates the +* multiplier coefficients for each haircell filter and +* initializes the filter states. It is called by function +* ``Meddis()'' in file ``model.c''. It returns a pointer to a +* structure containing all the relevant information for the +* computation of the haircell filters. +* +* The choice of the haircell type (medium or high-spontaneous +* rate fibre) is made at run time using options ``fibre_med''. +* +* The choice of the numerical implementation (Wave Digital +* Filtering or Forward Difference algorithm) is made at +* compile time using a macro substitution ``#define'' (see above). +* +* The choice of clamping the reservoirs' contents to non- +* negative values is made at compile time using a macro +* substitution ``#define'' (see above). +* +* DoHaircells() Realization of the bank of haircell filters. It computes the +* output for each haircell for a specified number of input +* points and keeps track of the filters' states. It returns +* the size of the output data in bytes (per point). +* +* CloseHaircells() Deletes all private data structures and arrays of the Haircell +* filters upon comletion of filtering. +********************************************************************************/ + + +/*************************** NewHaircells() *********************************/ + +Pointer NewHaircells( channels, samplerate, output_gain, coupling, fibreType ) + int channels ; + double samplerate, output_gain ; + double coupling ; + char *fibreType ; +{ + DeclareNew( HaircellBank *, bank ) ; + HaircellInfo *cell_info ; + CoeffType *cf ; + StateType *st ; + double ydt, ldt, rdt, xdt, gdt ; + double kdt_init, q_init, c_init, w_init, cell_scaling ; + double A_value, B_value, g_value ; + double r1, r2, r3, r4, r03 ; + double update_inputGain() ; + int channel ; + + /*** scaling of input data ***/ + cell_scaling = coupling ; + + /*** allocate fibre-dependent haircell parameters ***/ + if( strncmp( fibreType, "high", 3 ) == 0 ) { + A_value = A_high ; + B_value = B_high ; + g_value = g_high ; + } + else { + A_value = A_medium ; + B_value = B_medium ; + g_value = g_medium ; + } + +#ifdef _DEBUG_ + fprintf( stderr, "Meddis Haircell fibre type=%s\n", fibreType ) ; + fprintf( stderr, "A=%.2f B=%.2f g=%.2f\n", A_value, B_value, g_value ) ; + fprintf( stderr, "Input scaling factor=%.4f\n", cell_scaling ) ; +#endif + + + /*** scale model parameters for difference equations ***/ + ydt = y_value / samplerate ; + ldt = l_value / samplerate ; + rdt = r_value / samplerate ; + xdt = x_value / samplerate ; + gdt = g_value / samplerate ; + + + /*** calculate initial values for infinite history of silence ***/ + kdt_init = gdt * A_value / ( A_value + B_value ) ; + q_init = m_value / ( 1. + kdt_init / ydt * ( 1. - rdt / ( ldt + rdt ) ) ) ; + c_init = q_init * kdt_init / ( ldt + rdt ) ; + w_init = rdt / xdt * c_init ; + + + /*** allocate and store multiplier coefficients ***/ + bank->number_of_coeffs = NUMBER_OF_COEFFS ; + cf = bank->coeffs = NewArray( CoeffType, bank->number_of_coeffs, + "haircell_tl.c for coefficients" ) ; +#ifdef _WDF_ + output_gain = MEDDIS_SCALE * output_gain * h_value / ( 2. * r_value ) ; + cf[0] = m_value * y_value * output_gain ; /* M/Zy (normalized) */ + cf[1] = A_value / cell_scaling ; + cf[2] = ( A_value + B_value ) / cell_scaling ; + cf[3] = g_value ; + + /* Adaptor B0 */ + r1 = y_value ; /* ( 1/Zy ) */ + r2 = 2. * samplerate ; /* ( 1/Zcq ) */ + r3 = r03 = r1 + r2 ; + cf[5] = r1 / r3 ; /* gamma01 */ + cf[4] = r3 / cf[3] ; /* needed to update gamma11 */ + + /* Adaptor B1*/ + r1 = kdt_init * samplerate ; /* ( 1/Zk ) initial value only */ + r2 = r03 ; + r3 = r1 + r2 ; + cf[6] = 2. * r1 / r3 ; /* gamma11 updated at each sample */ + + /* Adaptor B2 */ + r1 = l_value ; /* ( 1/Zl ) */ + r2 = r_value ; /* ( 1/Zr ) */ + r3 = 2. * samplerate ; /* ( 1/Zcc ) */ + r4 = r1 + r2 + r3 ; + cf[7] = r1 / r4 ; /* gamma21 */ + cf[8] = r2 / r4 ; /* gamma22 */ + + /* Adaptor B3 */ + r1 = x_value ; /* ( 1/Zx ) */ + r2 = 2. * samplerate ; /* ( 1/Zcw ) */ + r3 = r1 + r2 ; + cf[9] = 0.5 * r1 / r3 ; /* 0.5 * gamma31 */ + +#else + output_gain = MEDDIS_SCALE * output_gain * h_value ; + cf[0] = m_value * output_gain ; + cf[1] = A_value / cell_scaling ; + cf[2] = ( A_value + B_value ) / cell_scaling ; + cf[3] = gdt ; + cf[4] = ydt ; + cf[5] = ldt ; + cf[6] = rdt ; + cf[7] = xdt ; +#endif + + + /*** loop through each haircell ***/ + bank->channels = channels ; + bank->cells = NewArray( HaircellInfo *, bank->channels, "haircell.c for cell states" ) ; + bank->inputGain_ratio = InputGain_ratio = NewArray( double, bank->channels, + "haircell_tl.c for cells input gain ratio" ) ; + bank->update_proc = update_inputGain ; + + for( channel = 0 ; channel < bank->channels ; channel++ ) { + bank->inputGain_ratio[ channel ] = 0.5 ; /* in absence of feedback, Gain = 1 */ + cell_info = bank->cells[ channel ] = New( HaircellInfo * ) ; + + /*** allocate and initialise states ***/ + cell_info->number_of_states = NUMBER_OF_STATES ; + st = cell_info->states = NewZeroedArray( StateType, cell_info->number_of_states, + "haircell_tl.c for states" ) ; + +#ifdef _WDF_ + st[1] = x_value * w_init * output_gain ; /* Ixn (normalized) */ + st[2] = 2. * q_init * samplerate * output_gain ; /* Ixn(t-T)-b02 (normalized) */ + st[3] = -2. * c_init * samplerate * output_gain ; /* b23 (normalized) */ + st[4] = -2. * w_init * samplerate * output_gain ; /* b33 (normalized) */ + cell_info->output_gain = output_gain ; /* obsolete */ + +#else + st[0] = kdt_init ; + st[1] = q_init * output_gain ; + st[2] = c_init * output_gain ; + st[3] = w_init * output_gain ; + cell_info->output_gain = output_gain ; /* obsolete */ +#endif + + } + + return( ( Pointer ) bank ) ; +} + + +/************************** DoHaircells() *********************************/ + +int DoHaircells( bank_ptr, bytes, output_ptr, end_output_ptr, input_ptr ) + Pointer *bank_ptr ; + ByteCount *bytes ; + DataType *output_ptr, *end_output_ptr, *input_ptr ; +{ + register HaircellBank *bank = ( HaircellBank * ) bank_ptr ; + register HaircellInfo *cell_info ; + register DataType *input, *output ; + register StateType *st ; + register StateType in, out ; + register CoeffType *cf ; + register double inputGain ; + register int channel = bank->channels ; + register int point, npoints = ( end_output_ptr - output_ptr ) / bank->channels ; +#ifdef _WDF_ + register StateType an, a0 ; +#else + register StateType replenish, eject, loss, reuptake, reprocess ; +#endif + + /*** for all channels ***/ + +#ifdef _DEBUG_ + fprintf( stderr, "DoHaircells() = %d points X %d channels\n", npoints, bank->channels ) ; +#endif + + cf = bank->coeffs ; + + while( channel-- ) { + + cell_info = bank->cells[ channel ] ; + st = cell_info->states ; + input = input_ptr + channel ; + output = output_ptr + channel ; + inputGain = bank->update_proc( bank->inputGain_ratio[ channel ] ) ; + cf[10] = cf[1] / inputGain ; /* (A/p) */ + cf[11] = cf[2] / inputGain ; /* (A+B)/p */ + + /*** for each channel ***/ + + point = npoints ; + while( point-- ) { + +#ifdef _WDF_ + /*** compute compressive permeability function k ***/ + in = ( StateType ) *input ; /* (BM vibration) */ + if( in > -cf[10] ) + st[0] = ( in + cf[10] ) / ( in + cf[11] ) ; /* kn(t) / g */ + else + st[0] = 0. ; + + /*** update gamma11 ***/ + cf[6] = st[0] / ( st[0] + cf[4] ) ; + cf[6] += cf[6] ; /* gamma11 */ + + /*** compute wdf ***/ + + /* Adaptor B0 */ + st[5] = cf[0] + st[1] + st[2] ; /* -b03=a12 (normalized) */ + + /* Adaptor B1 */ + st[6] = cf[6] * st[5] ; /* -b11=2Ikn (normalized) */ + st[7] = st[6] - st[5] ; /* b12=-a03 (normalized) */ + + /* Adaptor B0 */ + st[2] = cf[0] + st[1] - st[7] + cf[5] * ( st[7] - st[5] ) ; /* Ixn(t-T)-b02 (normalized) */ + + /* Adaptor B2 */ + an = st[6] - st[3] ; /* a24 (normalized) */ + a0 = an - st[3] ; /* a20 (normalized) */ + st[8] = cf[8] * a0 ; /* -b22=2Irn (normalized) */ + st[3] = cf[7] * a0 + st[8] - an ; /* b23 (normalized) */ + + /* Adaptor B3 */ + an = st[8] - st[4] ; /* a33 (normalized) */ + st[1] = cf[9] * ( an - st[4] ) ; /* Ixn(t)=-0.5*b31 (normalized) */ + st[4] = st[1] + st[1] - an ; /* b32 (normalized) */ + +#ifdef _CLAMPING_ + /*** clamping of reservoirs to non-negative quantities ***/ + /*** not available yet ***/ +#endif + + /*** output ***/ + out = st[8] ; + +#else + + /*** compute change quantities ***/ + replenish = cf[4] * ( cf[0] - st[1] ) ; /* y(M-q) */ + eject = st[0] * st[1] ; /* kq */ + loss = cf[5] * st[2] ; /* lc */ + reuptake = cf[6] * st[2] ; /* rc */ + reprocess = cf[7] * st[3] ; /* xw */ + + /*** compute compressive permeability function kn_1dt ***/ + in = ( StateType ) *input ; /* (BM vibration) */ + if( in > -cf[10] ) + st[0] = cf[3] * ( in + cf[10] ) / ( in + cf[11] ) ; /* kn_1(t)dt */ + else + st[0] = 0. ; + + /*** update reservoir quantities ***/ + st[1] += replenish - eject + reprocess ; /* q */ + st[2] += eject - loss - reuptake ; /* c */ + st[3] += reuptake - reprocess ; /* w */ + +#ifdef _CLAMPING_ + /*** clamping of reservoirs to non-negative quantities ***/ + if( st[1] < 0 ) + st[1] = 0 ; + if( st[2] < 0 ) + st[2] = 0 ; +#endif + + /*** output ***/ + out = st[2] ; + +#endif + + /*** output ***/ +#ifdef _OVERFLOW_ + if( out > _MaxOutput_ ) + fprintf( stderr, "Overflow error in Meddis Haircell output\%e\n", ( double ) out ) ; +#endif + + *output = ( DataType ) out ; + + output += bank->channels ; + input += bank->channels ; + } + + } + + return( npoints ) ; + } + + +/************************** CloseHaircells() *********************************/ + +void CloseHaircells( bank_ptr ) + Pointer bank_ptr ; +{ + HaircellBank *bank = ( HaircellBank * ) bank_ptr ; + HaircellInfo *cell_info ; + int channel ; + + for( channel = 0 ; channel < bank->channels ; channel++) { + Delete( bank->cells[ channel ]->states ) ; + Delete( bank->cells[ channel ] ) ; + } + + Delete( bank->cells ) ; + Delete( bank->coeffs ) ; + Delete( bank->inputGain_ratio ) ; + Delete( bank ) ; + + return ; +} + + +/********** low level functions ************/ + +double update_inputGain( ratio ) + double ratio ; +{ + return( pow( 10., ( ratio - 0.5 ) * InputGain_max / 20. ) ) ; +} +