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: 
tomwalters@0: 	formulae.c
tomwalters@0: 	==========
tomwalters@0: 
tomwalters@0:     Functions from : Moore and Glasberg, Formulae describing frequency
tomwalters@0:     selectivity as a function or frequency and level, and their use in
tomwalters@0:     calculating excitation paterns. Hearing reserch, 28 (1987) 209-225
tomwalters@0: 
tomwalters@0: */
tomwalters@0: 
tomwalters@0: #include <math.h>
tomwalters@0: #include <stdio.h>
tomwalters@0: #include "../glib/options.h"
tomwalters@0: 
tomwalters@0: #ifndef  _FORMULAE_H_
tomwalters@0: #include "formulae.h"
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: /* latest erb function circa. 1989 from BBG & CJM */
tomwalters@0: 
tomwalters@0: /* derived from erb = 24.7 * ( 4.37e-3 * f +  1 ) */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: static double limit = 24.7  ; /* limit = k1 */
tomwalters@0: static double Q     = 9.265 ; /* Q = 1. / ( k1 * k2 ) */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: void SetErbParameters( new_limit, new_Q )
tomwalters@0: double new_limit, new_Q ;
tomwalters@0: {
tomwalters@0:     limit = new_limit ;
tomwalters@0:     Q     = new_Q     ;
tomwalters@0: 
tomwalters@0:     return ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: double Erb( frequency )
tomwalters@0:      double frequency ;
tomwalters@0: {
tomwalters@0:     return ( limit + frequency / Q ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: /* integrate 1/erb(f) to get erb based scale */
tomwalters@0: 
tomwalters@0: double ErbScale( frequency )
tomwalters@0: 	  double frequency ;
tomwalters@0: {
tomwalters@0:     return ( log( 1. + frequency / Q / limit ) * Q ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: double FofErbScale( E )
tomwalters@0: 	     double E ;
tomwalters@0: {
tomwalters@0:   return ( ( exp( E / Q ) - 1 ) * Q * limit ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: double poly( coefts, x )
tomwalters@0: double *coefts, x ;
tomwalters@0: {
tomwalters@0:     double *cptr = coefts ;
tomwalters@0:     double power = 1 ;
tomwalters@0:     double value = 0 ;
tomwalters@0: 
tomwalters@0:     while( *cptr != 0. ) {
tomwalters@0: 	value += *cptr++ * power ;
tomwalters@0: 	power *= x ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     return ( value ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: double dBAudiogram( frequency )
tomwalters@0: 	     double frequency ;
tomwalters@0: {
tomwalters@0:     static double audiogram_polynomial_coeficients[] = {
tomwalters@0: 	  2661.8, -3690.1, 1917.4, -440.77, 37.706, 0. } ;
tomwalters@0: 
tomwalters@0:     return ( poly( audiogram_polynomial_coeficients, log10( frequency ) ) ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: double Audiogram( frequency )
tomwalters@0: 	   double frequency ;
tomwalters@0: {
tomwalters@0:     return ( pow( 10., -dBAudiogram( frequency ) / 20. ) ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /* old versions of above */
tomwalters@0: 
tomwalters@0: static double erb_a  = {  6.23e-6 },
tomwalters@0: 	      erb_b  = { 93.39e-3 },
tomwalters@0: 	      erb_c  = { 28.52    } ;
tomwalters@0: 
tomwalters@0: /* erb rate function coeficients */
tomwalters@0: 
tomwalters@0: static double erb_k1 = {  11.17   },
tomwalters@0: 	      erb_k2 = {   0.312  },
tomwalters@0: 	      erb_k3 = {  14.675  },
tomwalters@0: 	      erb_k4 = {  43.0    } ;
tomwalters@0: 
tomwalters@0: static double     b0 = {  2661.8  },
tomwalters@0: 		  b1 = {  3690.1  },
tomwalters@0: 		  b2 = {  1917.4  },
tomwalters@0: 		  b3 = {  440.77  },
tomwalters@0: 		  b4 = {  37.706  } ;
tomwalters@0: 
tomwalters@0: double OldErb( f )
tomwalters@0: double f ;
tomwalters@0: {
tomwalters@0:     return ( erb_a * f * f + erb_b * f + erb_c ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: double OldErbScale( f )
tomwalters@0: double f ;
tomwalters@0: {
tomwalters@0:     double fkHz ;
tomwalters@0: 
tomwalters@0:     fkHz = f / 1000. ;
tomwalters@0: 
tomwalters@0:     return ( erb_k1 * log( ( fkHz + erb_k2 ) / ( fkHz + erb_k3 ) ) + erb_k4 ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: double OldFofErbScale( E )
tomwalters@0: double E ;
tomwalters@0: {
tomwalters@0:     double tmp ;
tomwalters@0: 
tomwalters@0:     tmp = exp( ( E - erb_k4 ) / erb_k1 ) ;
tomwalters@0: 
tomwalters@0:     return ( ( erb_k2 - erb_k3 * tmp ) / ( tmp - 1.0 ) * 1000. ) ;
tomwalters@0: }
tomwalters@0: