comparison filter/formulae.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
26 formulae.c
27 ==========
28
29 Functions from : Moore and Glasberg, Formulae describing frequency
30 selectivity as a function or frequency and level, and their use in
31 calculating excitation paterns. Hearing reserch, 28 (1987) 209-225
32
33 */
34
35 #include <math.h>
36 #include <stdio.h>
37 #include "../glib/options.h"
38
39 #ifndef _FORMULAE_H_
40 #include "formulae.h"
41 #endif
42
43 /* latest erb function circa. 1989 from BBG & CJM */
44
45 /* derived from erb = 24.7 * ( 4.37e-3 * f + 1 ) */
46
47
48 static double limit = 24.7 ; /* limit = k1 */
49 static double Q = 9.265 ; /* Q = 1. / ( k1 * k2 ) */
50
51
52 void SetErbParameters( new_limit, new_Q )
53 double new_limit, new_Q ;
54 {
55 limit = new_limit ;
56 Q = new_Q ;
57
58 return ;
59 }
60
61
62 double Erb( frequency )
63 double frequency ;
64 {
65 return ( limit + frequency / Q ) ;
66 }
67
68 /* integrate 1/erb(f) to get erb based scale */
69
70 double ErbScale( frequency )
71 double frequency ;
72 {
73 return ( log( 1. + frequency / Q / limit ) * Q ) ;
74 }
75
76
77 double FofErbScale( E )
78 double E ;
79 {
80 return ( ( exp( E / Q ) - 1 ) * Q * limit ) ;
81 }
82
83
84 double poly( coefts, x )
85 double *coefts, x ;
86 {
87 double *cptr = coefts ;
88 double power = 1 ;
89 double value = 0 ;
90
91 while( *cptr != 0. ) {
92 value += *cptr++ * power ;
93 power *= x ;
94 }
95
96 return ( value ) ;
97 }
98
99 double dBAudiogram( frequency )
100 double frequency ;
101 {
102 static double audiogram_polynomial_coeficients[] = {
103 2661.8, -3690.1, 1917.4, -440.77, 37.706, 0. } ;
104
105 return ( poly( audiogram_polynomial_coeficients, log10( frequency ) ) ) ;
106 }
107
108 double Audiogram( frequency )
109 double frequency ;
110 {
111 return ( pow( 10., -dBAudiogram( frequency ) / 20. ) ) ;
112 }
113
114
115 /* old versions of above */
116
117 static double erb_a = { 6.23e-6 },
118 erb_b = { 93.39e-3 },
119 erb_c = { 28.52 } ;
120
121 /* erb rate function coeficients */
122
123 static double erb_k1 = { 11.17 },
124 erb_k2 = { 0.312 },
125 erb_k3 = { 14.675 },
126 erb_k4 = { 43.0 } ;
127
128 static double b0 = { 2661.8 },
129 b1 = { 3690.1 },
130 b2 = { 1917.4 },
131 b3 = { 440.77 },
132 b4 = { 37.706 } ;
133
134 double OldErb( f )
135 double f ;
136 {
137 return ( erb_a * f * f + erb_b * f + erb_c ) ;
138 }
139
140 double OldErbScale( f )
141 double f ;
142 {
143 double fkHz ;
144
145 fkHz = f / 1000. ;
146
147 return ( erb_k1 * log( ( fkHz + erb_k2 ) / ( fkHz + erb_k3 ) ) + erb_k4 ) ;
148 }
149
150 double OldFofErbScale( E )
151 double E ;
152 {
153 double tmp ;
154
155 tmp = exp( ( E - erb_k4 ) / erb_k1 ) ;
156
157 return ( ( erb_k2 - erb_k3 * tmp ) / ( tmp - 1.0 ) * 1000. ) ;
158 }
159