1 /*------------- Telecommunications & Signal Processing Lab -------------
2 McGill University
3
4 Routine:
5 void PQWOME (const double f[], double W2[], int N)
6 void PQenThresh (const double f[], double Et[], int N)
7 void PQexIndex (const double f[], double s[], int N);
8 void PQintNoise (const double f[], double EIN[], int N)
9 void PQtConst (double t100, double tmin, double Fs, const double f[],
10 double a[], double b[], int N)
11
12 Purpose:
13 Table generation routines PQevalAudio
14
15 Description:
16 These routines generate tables used in PQevalAudio.
17 PQWOME - Outer and middle ear weighting (squared) at given frequencies
18 PQenThresh - Energy thresholds at given frequencies
19 PQexIndex - Excitation index at given frequencies
20 PQintNoise - Internal noise energy at given frequencies
21 PQtConst - Generate coefficients for frequency dependent time constants
22
23 Parameters:
24 PQWOME
25 -> const double f[]
26 Array of N frequency values
27 <- double W2[]
28 Array of N squared weight values
29 -> int N
30 Number of array elements
31 PQenThresh
32 -> const double f[]
33 Array of N frequency values
34 <- double Et[]
35 Array of N energy theshold values
36 -> int N
37 Number of array elements
38 PQexIndex
39 -> const double f[]
40 Array of N frequency values
41 <- double s[]
42 Array of N excitation index (energy) values
43 -> int N
44 Number of array elements
45 PQintNoise
46 -> const double f[]
47 Array of N frequency values
48 <- double EIN[]
49 Array of N internal noise energies
50 -> int N
51 Number of array elements
52 PQtConst
53 -> double t100
54 Time constant at 100 Hz (s)
55 -> double tmin
56 Minimum time constant (s)
57 -> double Fs
58 Sampling frequency (Hz)
59 -> const double f[]
60 Frequency array
61 <- double a[]
62 Difference equation coefficients (N values)
63 <- double b[]
64 Difference equation coefficients, b[k] = (1 - a[k]). b can be NULL.
65 -> int N
66 Number of frequency values and coefficient values
67
68 Author / revision:
69 P. Kabal Copyright (C) 2003
70 $Revision: 1.9 $ $Date: 2003/05/13 01:12:06 $
71
72 ----------------------------------------------------------------------*/
73
74 #include <math.h>
75
76 #include "PQevalAudio.h"
77
78 #define DBPI(x) (pow (10., x / 10.)) /* Inverse dB (power) */
79 #define SQRV(x) ((x) * (x))
80 #define PI 3.14159265358979323846
81
82 /*-------------------------------*/
83 /* Outer & middle ear weighting */
84
85 void
PQWOME(const double f[],double W2[],int N)86 PQWOME (const double f[], double W2[], int N)
87
88 {
89 int k;
90 double fkHz, AdB;
91
92 for (k = 0; k < N; ++k) {
93 if (f[k] != 0) {
94 fkHz = f[k] / 1000.;
95 AdB = -2.184 * pow (fkHz, -0.8) + 6.5 * exp(-0.6 * SQRV (fkHz - 3.3))
96 -0.001 * pow (fkHz, 3.6);
97 W2[k] = DBPI (AdB);
98 }
99 else {
100 W2[k] = 0;
101 }
102 }
103
104 return;
105 }
106
107 /* Energy threshold */
108
109 void
PQenThresh(const double f[],double Et[],int N)110 PQenThresh (const double f[], double Et[], int N)
111
112 {
113 int m;
114 double EtdB;
115
116 for (m = 0; m < N; ++m) {
117 EtdB = 3.64 * pow (f[m] / 1000., -0.8);
118 Et[m] = DBPI (EtdB);
119 }
120
121 return;
122 }
123
124 /* Excitation index */
125
126 void
PQexIndex(const double f[],double s[],int N)127 PQexIndex (const double f[], double s[], int N)
128
129 {
130 int m;
131 double sdB;
132
133 for (m = 0; m < N; ++m) {
134 sdB = -2 - 2.05 * atan (f[m]/4000.) - 0.75 * atan (SQRV(f[m]/1600.));
135 s[m] = DBPI (sdB);
136 }
137
138 return;
139 }
140
141 /* Internal noise */
142
143 void
PQintNoise(const double f[],double EIN[],int N)144 PQintNoise (const double f[], double EIN[], int N)
145
146 {
147 int m;
148 double INdB;
149
150 for (m = 0; m < N; ++m) {
151 INdB = 1.456 * pow (f[m] / 1000., -0.8);
152 EIN[m] = DBPI (INdB);
153 }
154
155 return;
156 }
157
158 /* Frequency dependent time constants */
159
160 void
PQtConst(double t100,double tmin,double Fs,const double f[],double a[],double b[],int N)161 PQtConst (double t100, double tmin, double Fs, const double f[],
162 double a[], double b[], int N)
163
164 {
165 int m;
166 double t;
167
168 for (m = 0; m < N; ++m) {
169 t = tmin + (100 / f[m]) * (t100 - tmin);
170 a[m] = exp (-1. / (Fs * t));
171 if (b != NULL)
172 b[m] = 1. - a[m];
173 }
174
175 return;
176 }
177