diff filter/formulae.c @ 0:5242703e91d3 tip

Initial checkin for AIM92 aimR8.2 (last updated May 1997).

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