Setting i to 0 once should do.
[ffmpeg.git] / libavcodec / fft.c
1 /*
2  * FFT/IFFT transforms
3  * Copyright (c) 2002 Fabrice Bellard.
4  *
5  * This file is part of FFmpeg.
6  *
7  * FFmpeg is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public
9  * License as published by the Free Software Foundation; either
10  * version 2.1 of the License, or (at your option) any later version.
11  *
12  * FFmpeg is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15  * Lesser General Public License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with FFmpeg; if not, write to the Free Software
19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20  */
21
22 /**
23  * @file fft.c
24  * FFT/IFFT transforms.
25  */
26
27 #include "dsputil.h"
28
29 /**
30  * The size of the FFT is 2^nbits. If inverse is TRUE, inverse FFT is
31  * done
32  */
33 int ff_fft_init(FFTContext *s, int nbits, int inverse)
34 {
35     int i, j, m, n;
36     float alpha, c1, s1, s2;
37     int shuffle = 0;
38     int av_unused has_vectors;
39
40     s->nbits = nbits;
41     n = 1 << nbits;
42
43     s->exptab = av_malloc((n / 2) * sizeof(FFTComplex));
44     if (!s->exptab)
45         goto fail;
46     s->revtab = av_malloc(n * sizeof(uint16_t));
47     if (!s->revtab)
48         goto fail;
49     s->inverse = inverse;
50
51     s2 = inverse ? 1.0 : -1.0;
52
53     for(i=0;i<(n/2);i++) {
54         alpha = 2 * M_PI * (float)i / (float)n;
55         c1 = cos(alpha);
56         s1 = sin(alpha) * s2;
57         s->exptab[i].re = c1;
58         s->exptab[i].im = s1;
59     }
60     s->fft_calc = ff_fft_calc_c;
61     s->imdct_calc = ff_imdct_calc;
62     s->exptab1 = NULL;
63
64 #ifdef HAVE_MMX
65     has_vectors = mm_support();
66     shuffle = 1;
67     if (has_vectors & MM_3DNOWEXT) {
68         /* 3DNowEx for K7/K8 */
69         s->imdct_calc = ff_imdct_calc_3dn2;
70         s->fft_calc = ff_fft_calc_3dn2;
71     } else if (has_vectors & MM_3DNOW) {
72         /* 3DNow! for K6-2/3 */
73         s->fft_calc = ff_fft_calc_3dn;
74     } else if (has_vectors & MM_SSE) {
75         /* SSE for P3/P4 */
76         s->imdct_calc = ff_imdct_calc_sse;
77         s->fft_calc = ff_fft_calc_sse;
78     } else {
79         shuffle = 0;
80     }
81 #elif defined HAVE_ALTIVEC && !defined ALTIVEC_USE_REFERENCE_C_CODE
82     has_vectors = mm_support();
83     if (has_vectors & MM_ALTIVEC) {
84         s->fft_calc = ff_fft_calc_altivec;
85         shuffle = 1;
86     }
87 #endif
88
89     /* compute constant table for HAVE_SSE version */
90     if (shuffle) {
91         int np, nblocks, np2, l;
92         FFTComplex *q;
93
94         np = 1 << nbits;
95         nblocks = np >> 3;
96         np2 = np >> 1;
97         s->exptab1 = av_malloc(np * 2 * sizeof(FFTComplex));
98         if (!s->exptab1)
99             goto fail;
100         q = s->exptab1;
101         do {
102             for(l = 0; l < np2; l += 2 * nblocks) {
103                 *q++ = s->exptab[l];
104                 *q++ = s->exptab[l + nblocks];
105
106                 q->re = -s->exptab[l].im;
107                 q->im = s->exptab[l].re;
108                 q++;
109                 q->re = -s->exptab[l + nblocks].im;
110                 q->im = s->exptab[l + nblocks].re;
111                 q++;
112             }
113             nblocks = nblocks >> 1;
114         } while (nblocks != 0);
115         av_freep(&s->exptab);
116     }
117
118     /* compute bit reverse table */
119
120     for(i=0;i<n;i++) {
121         m=0;
122         for(j=0;j<nbits;j++) {
123             m |= ((i >> j) & 1) << (nbits-j-1);
124         }
125         s->revtab[i]=m;
126     }
127     return 0;
128  fail:
129     av_freep(&s->revtab);
130     av_freep(&s->exptab);
131     av_freep(&s->exptab1);
132     return -1;
133 }
134
135 /* butter fly op */
136 #define BF(pre, pim, qre, qim, pre1, pim1, qre1, qim1) \
137 {\
138   FFTSample ax, ay, bx, by;\
139   bx=pre1;\
140   by=pim1;\
141   ax=qre1;\
142   ay=qim1;\
143   pre = (bx + ax);\
144   pim = (by + ay);\
145   qre = (bx - ax);\
146   qim = (by - ay);\
147 }
148
149 #define MUL16(a,b) ((a) * (b))
150
151 #define CMUL(pre, pim, are, aim, bre, bim) \
152 {\
153    pre = (MUL16(are, bre) - MUL16(aim, bim));\
154    pim = (MUL16(are, bim) + MUL16(bre, aim));\
155 }
156
157 /**
158  * Do a complex FFT with the parameters defined in ff_fft_init(). The
159  * input data must be permuted before with s->revtab table. No
160  * 1.0/sqrt(n) normalization is done.
161  */
162 void ff_fft_calc_c(FFTContext *s, FFTComplex *z)
163 {
164     int ln = s->nbits;
165     int j, np, np2;
166     int nblocks, nloops;
167     register FFTComplex *p, *q;
168     FFTComplex *exptab = s->exptab;
169     int l;
170     FFTSample tmp_re, tmp_im;
171
172     np = 1 << ln;
173
174     /* pass 0 */
175
176     p=&z[0];
177     j=(np >> 1);
178     do {
179         BF(p[0].re, p[0].im, p[1].re, p[1].im,
180            p[0].re, p[0].im, p[1].re, p[1].im);
181         p+=2;
182     } while (--j != 0);
183
184     /* pass 1 */
185
186
187     p=&z[0];
188     j=np >> 2;
189     if (s->inverse) {
190         do {
191             BF(p[0].re, p[0].im, p[2].re, p[2].im,
192                p[0].re, p[0].im, p[2].re, p[2].im);
193             BF(p[1].re, p[1].im, p[3].re, p[3].im,
194                p[1].re, p[1].im, -p[3].im, p[3].re);
195             p+=4;
196         } while (--j != 0);
197     } else {
198         do {
199             BF(p[0].re, p[0].im, p[2].re, p[2].im,
200                p[0].re, p[0].im, p[2].re, p[2].im);
201             BF(p[1].re, p[1].im, p[3].re, p[3].im,
202                p[1].re, p[1].im, p[3].im, -p[3].re);
203             p+=4;
204         } while (--j != 0);
205     }
206     /* pass 2 .. ln-1 */
207
208     nblocks = np >> 3;
209     nloops = 1 << 2;
210     np2 = np >> 1;
211     do {
212         p = z;
213         q = z + nloops;
214         for (j = 0; j < nblocks; ++j) {
215             BF(p->re, p->im, q->re, q->im,
216                p->re, p->im, q->re, q->im);
217
218             p++;
219             q++;
220             for(l = nblocks; l < np2; l += nblocks) {
221                 CMUL(tmp_re, tmp_im, exptab[l].re, exptab[l].im, q->re, q->im);
222                 BF(p->re, p->im, q->re, q->im,
223                    p->re, p->im, tmp_re, tmp_im);
224                 p++;
225                 q++;
226             }
227
228             p += nloops;
229             q += nloops;
230         }
231         nblocks = nblocks >> 1;
232         nloops = nloops << 1;
233     } while (nblocks != 0);
234 }
235
236 /**
237  * Do the permutation needed BEFORE calling ff_fft_calc()
238  */
239 void ff_fft_permute(FFTContext *s, FFTComplex *z)
240 {
241     int j, k, np;
242     FFTComplex tmp;
243     const uint16_t *revtab = s->revtab;
244
245     /* reverse */
246     np = 1 << s->nbits;
247     for(j=0;j<np;j++) {
248         k = revtab[j];
249         if (k < j) {
250             tmp = z[k];
251             z[k] = z[j];
252             z[j] = tmp;
253         }
254     }
255 }
256
257 void ff_fft_end(FFTContext *s)
258 {
259     av_freep(&s->revtab);
260     av_freep(&s->exptab);
261     av_freep(&s->exptab1);
262 }
263