7e95317ef53ede6ac4fb0fe6e65a03869934d002
[asterisk/asterisk.git] / codecs / mp3 / src / mdct.c
1 /*____________________________________________________________________________
2         
3         FreeAmp - The Free MP3 Player
4
5         MP3 Decoder originally Copyright (C) 1995-1997 Xing Technology
6         Corp.  http://www.xingtech.com
7
8         Portions Copyright (C) 1998-1999 EMusic.com
9
10         This program is free software; you can redistribute it and/or modify
11         it under the terms of the GNU General Public License as published by
12         the Free Software Foundation; either version 2 of the License, or
13         (at your option) any later version.
14
15         This program is distributed in the hope that it will be useful,
16         but WITHOUT ANY WARRANTY; without even the implied warranty of
17         MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18         GNU General Public License for more details.
19
20         You should have received a copy of the GNU General Public License
21         along with this program; if not, write to the Free Software
22         Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
23         
24         $Id$
25 ____________________________________________________________________________*/
26
27 /****  mdct.c  ***************************************************
28
29 Layer III 
30
31   cos transform for n=18, n=6
32
33 computes  c[k] =  Sum( cos((pi/4*n)*(2*k+1)*(2*p+1))*f[p] )
34                 k = 0, ...n-1,  p = 0...n-1
35
36
37 inplace ok.
38
39 ******************************************************************/
40
41 #include <stdlib.h>
42 #include <stdio.h>
43 #include <float.h>
44 #include <math.h>
45
46 #ifdef ASM_X86
47 extern void imdct18_asm(float f[18]);
48 extern void imdct6_3_asm(float f[]);
49 #endif /* ASM_X86 */
50
51 /*------ 18 point xform -------*/
52 float mdct18w[18];
53 float mdct18w2[9];
54 float coef[9][4];
55
56 float mdct6_3v[6];
57 float mdct6_3v2[3];
58 float coef87;
59
60 typedef struct
61 {
62    float *w;
63    float *w2;
64    void *coef;
65 }
66 IMDCT_INIT_BLOCK;
67
68 static IMDCT_INIT_BLOCK imdct_info_18 =
69 {mdct18w, mdct18w2, coef};
70 static IMDCT_INIT_BLOCK imdct_info_6 =
71 {mdct6_3v, mdct6_3v2, &coef87};
72
73 /*====================================================================*/
74 IMDCT_INIT_BLOCK *imdct_init_addr_18()
75 {
76    return &imdct_info_18;
77 }
78 IMDCT_INIT_BLOCK *imdct_init_addr_6()
79 {
80    return &imdct_info_6;
81 }
82 /*--------------------------------------------------------------------*/
83 void imdct18(float f[18])       /* 18 point */
84 {
85 #ifdef ASM_X86
86    imdct18_asm(f);
87 #else
88    int p;
89    float a[9], b[9];
90    float ap, bp, a8p, b8p;
91    float g1, g2;
92
93
94    for (p = 0; p < 4; p++)
95    {
96       g1 = mdct18w[p] * f[p];
97       g2 = mdct18w[17 - p] * f[17 - p];
98       ap = g1 + g2;             // a[p]
99
100       bp = mdct18w2[p] * (g1 - g2);     // b[p]
101
102       g1 = mdct18w[8 - p] * f[8 - p];
103       g2 = mdct18w[9 + p] * f[9 + p];
104       a8p = g1 + g2;            // a[8-p]
105
106       b8p = mdct18w2[8 - p] * (g1 - g2);        // b[8-p]
107
108       a[p] = ap + a8p;
109       a[5 + p] = ap - a8p;
110       b[p] = bp + b8p;
111       b[5 + p] = bp - b8p;
112    }
113    g1 = mdct18w[p] * f[p];
114    g2 = mdct18w[17 - p] * f[17 - p];
115    a[p] = g1 + g2;
116    b[p] = mdct18w2[p] * (g1 - g2);
117
118
119    f[0] = 0.5f * (a[0] + a[1] + a[2] + a[3] + a[4]);
120    f[1] = 0.5f * (b[0] + b[1] + b[2] + b[3] + b[4]);
121
122    f[2] = coef[1][0] * a[5] + coef[1][1] * a[6] + coef[1][2] * a[7]
123       + coef[1][3] * a[8];
124    f[3] = coef[1][0] * b[5] + coef[1][1] * b[6] + coef[1][2] * b[7]
125       + coef[1][3] * b[8] - f[1];
126    f[1] = f[1] - f[0];
127    f[2] = f[2] - f[1];
128
129    f[4] = coef[2][0] * a[0] + coef[2][1] * a[1] + coef[2][2] * a[2]
130       + coef[2][3] * a[3] - a[4];
131    f[5] = coef[2][0] * b[0] + coef[2][1] * b[1] + coef[2][2] * b[2]
132       + coef[2][3] * b[3] - b[4] - f[3];
133    f[3] = f[3] - f[2];
134    f[4] = f[4] - f[3];
135
136    f[6] = coef[3][0] * (a[5] - a[7] - a[8]);
137    f[7] = coef[3][0] * (b[5] - b[7] - b[8]) - f[5];
138    f[5] = f[5] - f[4];
139    f[6] = f[6] - f[5];
140
141    f[8] = coef[4][0] * a[0] + coef[4][1] * a[1] + coef[4][2] * a[2]
142       + coef[4][3] * a[3] + a[4];
143    f[9] = coef[4][0] * b[0] + coef[4][1] * b[1] + coef[4][2] * b[2]
144       + coef[4][3] * b[3] + b[4] - f[7];
145    f[7] = f[7] - f[6];
146    f[8] = f[8] - f[7];
147
148    f[10] = coef[5][0] * a[5] + coef[5][1] * a[6] + coef[5][2] * a[7]
149       + coef[5][3] * a[8];
150    f[11] = coef[5][0] * b[5] + coef[5][1] * b[6] + coef[5][2] * b[7]
151       + coef[5][3] * b[8] - f[9];
152    f[9] = f[9] - f[8];
153    f[10] = f[10] - f[9];
154
155    f[12] = 0.5f * (a[0] + a[2] + a[3]) - a[1] - a[4];
156    f[13] = 0.5f * (b[0] + b[2] + b[3]) - b[1] - b[4] - f[11];
157    f[11] = f[11] - f[10];
158    f[12] = f[12] - f[11];
159
160    f[14] = coef[7][0] * a[5] + coef[7][1] * a[6] + coef[7][2] * a[7]
161       + coef[7][3] * a[8];
162    f[15] = coef[7][0] * b[5] + coef[7][1] * b[6] + coef[7][2] * b[7]
163       + coef[7][3] * b[8] - f[13];
164    f[13] = f[13] - f[12];
165    f[14] = f[14] - f[13];
166
167    f[16] = coef[8][0] * a[0] + coef[8][1] * a[1] + coef[8][2] * a[2]
168       + coef[8][3] * a[3] + a[4];
169    f[17] = coef[8][0] * b[0] + coef[8][1] * b[1] + coef[8][2] * b[2]
170       + coef[8][3] * b[3] + b[4] - f[15];
171    f[15] = f[15] - f[14];
172    f[16] = f[16] - f[15];
173    f[17] = f[17] - f[16];
174
175
176    return;
177 #endif
178 }
179 /*--------------------------------------------------------------------*/
180 /* does 3, 6 pt dct.  changes order from f[i][window] c[window][i] */
181 void imdct6_3(float f[])        /* 6 point */
182 {
183 #ifdef ASM_X86
184    imdct6_3_asm(f);
185 #else
186    int w;
187    float buf[18];
188    float *a, *c;                // b[i] = a[3+i]
189
190    float g1, g2;
191    float a02, b02;
192
193    c = f;
194    a = buf;
195    for (w = 0; w < 3; w++)
196    {
197       g1 = mdct6_3v[0] * f[3 * 0];
198       g2 = mdct6_3v[5] * f[3 * 5];
199       a[0] = g1 + g2;
200       a[3 + 0] = mdct6_3v2[0] * (g1 - g2);
201
202       g1 = mdct6_3v[1] * f[3 * 1];
203       g2 = mdct6_3v[4] * f[3 * 4];
204       a[1] = g1 + g2;
205       a[3 + 1] = mdct6_3v2[1] * (g1 - g2);
206
207       g1 = mdct6_3v[2] * f[3 * 2];
208       g2 = mdct6_3v[3] * f[3 * 3];
209       a[2] = g1 + g2;
210       a[3 + 2] = mdct6_3v2[2] * (g1 - g2);
211
212       a += 6;
213       f++;
214    }
215
216    a = buf;
217    for (w = 0; w < 3; w++)
218    {
219       a02 = (a[0] + a[2]);
220       b02 = (a[3 + 0] + a[3 + 2]);
221       c[0] = a02 + a[1];
222       c[1] = b02 + a[3 + 1];
223       c[2] = coef87 * (a[0] - a[2]);
224       c[3] = coef87 * (a[3 + 0] - a[3 + 2]) - c[1];
225       c[1] = c[1] - c[0];
226       c[2] = c[2] - c[1];
227       c[4] = a02 - a[1] - a[1];
228       c[5] = b02 - a[3 + 1] - a[3 + 1] - c[3];
229       c[3] = c[3] - c[2];
230       c[4] = c[4] - c[3];
231       c[5] = c[5] - c[4];
232       a += 6;
233       c += 6;
234    }
235
236    return;
237 #endif
238 }
239 /*--------------------------------------------------------------------*/