1 /*
2 * Musepack audio compression
3 * Copyright (c) 2005-2009, The Musepack Development Team
4 * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf
5 *
6 * This library is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 of the License, or (at your option) any later version.
10 *
11 * This library is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with this library; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 */
20
21 #include <string.h>
22
23 #include "libmpcpsy.h"
24 #include <mpc/mpcmath.h>
25
26 #define CX0 -1.
27 #define CX1 0.5
28
29 #define SX1 -1.
30 #define SX2 (2./9/ 1)
31 #define SX3 (2./9/ 4)
32 #define SX4 (2./9/ 10)
33 #define SX5 (2./9/ 20)
34 #define SX6 (2./9/ 35)
35 #define SX7 (2./9/ 56)
36 #define SX8 (2./9/ 84)
37 #define SX9 (2./9/120)
38 #define SX10 (2./9/165)
39
40
41 #ifdef EXTRA_DECONV
42 # define DECONV \
43 { \
44 tmp = (CX0*aix[0] + CX1*aix[2]) * (1./(CX0*CX0+CX1*CX1)); \
45 aix[ 0] -= CX0*tmp; \
46 aix[ 2] -= CX1*tmp; \
47 tmp = (SX1*aix[3] + SX2*aix[5] + SX3*aix[7] + SX4*aix[9] + SX5*aix[11]) * (1./(SX1*SX1+SX2*SX2+SX3*SX3+SX4*SX4+SX5*SX5)); \
48 aix[ 3] -= SX1*tmp; \
49 aix[ 5] -= SX2*tmp; \
50 aix[ 7] -= SX3*tmp; \
51 aix[ 9] -= SX4*tmp; \
52 aix[11] -= SX5*tmp; \
53 }
54 #elif 0
55 # define DECONV \
56 { \
57 float A[20]; \
58 int i; \
59 memcpy (A, aix, 20*sizeof(aix)); \
60 tmp = (CX0*aix[0] + CX1*aix[2]) * (1./(CX0*CX0+CX1*CX1)); \
61 aix[ 0] -= CX0*tmp; \
62 aix[ 2] -= CX1*tmp; \
63 tmp = (SX1*aix[3] + SX2*aix[5] + SX3*aix[7] + SX4*aix[9] + SX5*aix[11]) * (1./(SX1*SX1+SX2*SX2+SX3*SX3+SX4*SX4+SX5*SX5)); \
64 aix[ 3] -= SX1*tmp; \
65 aix[ 5] -= SX2*tmp; \
66 aix[ 7] -= SX3*tmp; \
67 aix[ 9] -= SX4*tmp; \
68 aix[11] -= SX5*tmp; \
69 for ( i=0; i<10; i++) \
70 printf ("%u%9.0f%7.0f%9.0f%7.0f\n",i, A[i+i], A[i+i+1], aix[i+i], aix[i+i+1] ); \
71 }
72 #else
73 # define DECONV
74 #endif
75
76
77 /* V A R I A B L E S */
78 static int ip [4096]; // bitinverse for maximum 2048 FFT
79 static float w [4096]; // butterfly-coefficient for maximum 2048 FFT
80 static float a [4096]; // holds real input for FFT
81 static float Hann_256 [ 256];
82 static float Hann_1024 [1024];
83 static float Hann_1600 [1600];
84
85 void Generate_FFT_Tables ( const int, int*, float* );
86 void rdft ( const int, float*, int*, float* );
87
88
89 //////////////////////////////
90 //
91 // BesselI0 -- Regular Modified Cylindrical Bessel Function (Bessel I).
92 //
93
94 static double
Bessel_I_0(double x)95 Bessel_I_0 ( double x )
96 {
97 double denominator;
98 double numerator;
99 double z;
100
101 if (x == 0.)
102 return 1.;
103
104 z = x * x;
105 numerator = z* (z* (z* (z* (z* (z* (z* (z* (z* (z* (z* (z* (z* (z*
106 0.210580722890567e-22 + 0.380715242345326e-19 ) +
107 0.479440257548300e-16) + 0.435125971262668e-13 ) +
108 0.300931127112960e-10) + 0.160224679395361e-07 ) +
109 0.654858370096785e-05) + 0.202591084143397e-02 ) +
110 0.463076284721000e+00) + 0.754337328948189e+02 ) +
111 0.830792541809429e+04) + 0.571661130563785e+06 ) +
112 0.216415572361227e+08) + 0.356644482244025e+09 ) +
113 0.144048298227235e+10;
114
115 denominator = z* (z* (z - 0.307646912682801e+04) + 0.347626332405882e+07) - 0.144048298227235e+10;
116
117 return - numerator / denominator;
118 }
119
120 static double
residual(double x)121 residual ( double x )
122 {
123 return sqrt ( 1. - x*x );
124 }
125
126 //////////////////////////////
127 //
128 // KBDWindow -- Kaiser Bessel Derived Window
129 // fills the input window array with size samples of the
130 // KBD window with the given tuning parameter alpha.
131 //
132
133
134 static void
KBDWindow(float * window,unsigned int size,float alpha)135 KBDWindow ( float* window, unsigned int size, float alpha )
136 {
137 double sumvalue = 0.;
138 double scale;
139 int i;
140
141 scale = 0.25 / sqrt (size);
142 for ( i = 0; i < (int)size/2; i++ )
143 window [i] = sumvalue += Bessel_I_0 ( M_PI * alpha * residual (4.*i/size - 1.) );
144
145 // need to add one more value to the nomalization factor at size/2:
146 sumvalue += Bessel_I_0 ( M_PI * alpha * residual (4.*(size/2)/size-1.) );
147
148 // normalize the window and fill in the righthand side of the window:
149 for ( i = 0; i < (int)size/2; i++ )
150 window [size-1-i] = window [i] = /*sqrt*/ ( window [i] / sumvalue ) * scale;
151 }
152
153 static void
CosWindow(float * window,unsigned int size)154 CosWindow ( float* window, unsigned int size )
155 {
156 double x;
157 double scale;
158 int i;
159
160 scale = 0.25 / sqrt (size);
161 for ( i = 0; i < (int)size/2; i++ ) {
162 x = cos ( (i+0.5) * (M_PI / size) );
163 window [size/2-1-i] = window [size/2+i] = scale * x * x;
164 }
165 }
166
167 static void
Window(float * window,unsigned int size,float alpha)168 Window ( float* window, unsigned int size, float alpha )
169 {
170 if ( alpha < 0. )
171 CosWindow ( window, size ) ;
172 else
173 KBDWindow ( window, size, alpha );
174 }
175
176 /* F U N C T I O N S */
177 // generates FFT lookup-tables
178 void
Init_FFT(PsyModel * m)179 Init_FFT ( PsyModel* m )
180 {
181 int n;
182 double x;
183 double scale;
184
185 // normalized hann functions
186 Window ( Hann_256 , 256, m->KBD1 );
187 Window ( Hann_1024, 1024, m->KBD2 );
188 scale = 0.25 / sqrt (2048.);
189 for ( n = 0; n < 800; n++ )
190 x = cos ((n+0.5) * (M_PI/1600)), Hann_1600 [799-n] = Hann_1600 [800+n] = (float)(x * x * scale);
191
192 Generate_FFT_Tables ( 2048, ip, w );
193 }
194
195 // input : Signal *x
196 // output: energy spectrum *erg
197 void
PowSpec256(const float * x,float * erg)198 PowSpec256 ( const float* x, float* erg )
199 {
200 int i;
201 // windowing
202 for( i = 0; i < 256; ++i )
203 a[i] = x[i] * Hann_256[i];
204
205 rdft(256, a, ip, w); // perform FFT
206
207 // calculate power
208 for( i = 0; i < 128; ++i)
209 erg[i] = a[i*2] * a[i*2] + a[i*2+1] * a[i*2+1];
210 }
211
212 // input : Signal *x
213 // output: energy spectrum *erg
214 void
PowSpec1024(const float * x,float * erg)215 PowSpec1024 ( const float* x, float* erg )
216 {
217 int i;
218 // windowing
219 for( i = 0; i < 1024; ++i )
220 a[i] = x[i] * Hann_1024[i];
221
222 rdft(1024, a, ip, w); // perform FFT
223
224 // calculate power
225 for( i = 0; i < 512; ++i)
226 erg[i] = a[i*2] * a[i*2] + a[i*2+1] * a[i*2+1];
227 }
228
229 // input : Signal *x
230 // output: energy spectrum *erg
231 void
PowSpec2048(const float * x,float * erg)232 PowSpec2048 ( const float* x, float* erg )
233 {
234 int i;
235 // windowing (only 1600 samples available -> centered in 2048!)
236 memset(a, 0, 224 * sizeof *a);
237 for( i = 0; i < 1600; ++i )
238 a[i+224] = x[i] * Hann_1600[i];
239 memset(a + 1824, 0, 224 * sizeof *a);
240
241 rdft(2048, a, ip, w); // perform FFT
242
243 // calculate power
244 for( i = 0; i < 1024; ++i)
245 erg[i] = a[i*2] * a[i*2] + a[i*2+1] * a[i*2+1];
246 }
247
248 // input : Signal *x
249 // output: energy spectrum *erg and phase spectrum *phs
250 void
PolarSpec1024(const float * x,float * erg,float * phs)251 PolarSpec1024 ( const float* x, float* erg, float* phs )
252 {
253 int i;
254 for( i = 0; i < 1024; i++ )
255 a[i] = x[i] * Hann_1024[i];
256
257 rdft( 1024, a, ip, w); // perform FFT
258
259 // calculate power and phase
260 for( i = 0; i < 512; ++i )
261 {
262 erg[i] = a[i*2] * a[i*2] + a[i*2+1] * a[i*2+1];
263 phs[i] = ATAN2F( a[i*2+1], a[i*2] );
264 }
265 }
266
267 // input : logarithmized energy spectrum *cep
268 // output: Cepstrum *cep (in-place)
269 void
Cepstrum2048(float * cep,const int MaxLine)270 Cepstrum2048 ( float* cep, const int MaxLine )
271 {
272 int i, j;
273 // generate real, even spectrum (symmetric around 1024, cep[2048-i] = cep[i])
274 for( i = 0, j = 1024; i < 1024; ++i, --j )
275 cep[1024 + j] = cep[i];
276
277 rdft(2048, cep, ip, w);
278
279 // only real part as outcome (all even indexes of cep[])
280 for( i = 0; i < MaxLine + 1; ++i )
281 cep[i] = cep[i*2] * (float) (0.9888 / 2048);
282 }
283