4abdd4f7403470b1a22e3dd395170ad8c4b99667
[ffmpeg.git] / libavcodec / ppc / fft_altivec.c
1 /*
2  * FFT/IFFT transforms
3  * AltiVec-enabled
4  * Copyright (c) 2003 Romain Dolbeau <romain@dolbeau.org>
5  * Based on code Copyright (c) 2002 Fabrice Bellard.
6  *
7  * This library 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 of the License, or (at your option) any later version.
11  *
12  * This library 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 this library; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
20  */
21 #include "../dsputil.h"
22
23 #include "dsputil_altivec.h"
24
25 /**
26  * Do a complex FFT with the parameters defined in fft_init(). The
27  * input data must be permuted before with s->revtab table. No
28  * 1.0/sqrt(n) normalization is done.
29  * AltiVec-enabled
30  * This code assumes that the 'z' pointer is 16 bytes-aligned
31  * It also assumes all FFTComplex are 8 bytes-aligned pair of float
32  * The code is exactly the same as the SSE version, except
33  * that successive MUL + ADD/SUB have been merged into
34  * fused multiply-add ('vec_madd' in altivec)
35  */
36 void fft_calc_altivec(FFTContext *s, FFTComplex *z)
37 {
38     register const vector float vczero = (const vector float)(0.);
39     
40     int ln = s->nbits;
41     int j, np, np2;
42     int nblocks, nloops;
43     register FFTComplex *p, *q;
44     FFTComplex *cptr, *cptr1;
45     int k;
46
47     np = 1 << ln;
48
49     {
50         vector float *r, a, b, a1, c1, c2;
51
52         r = (vector float *)&z[0];
53
54         c1 = vcii(p,p,n,n);
55         
56         if (s->inverse)
57             {
58                 c2 = vcii(p,p,n,p);
59             }
60         else
61             {
62                 c2 = vcii(p,p,p,n);
63             }
64         
65         j = (np >> 2);
66         do {
67             a = vec_ld(0, r);
68             a1 = vec_ld(sizeof(vector float), r);
69             
70             b = vec_perm(a,a,vcprmle(1,0,3,2));
71             a = vec_madd(a,c1,b);
72             /* do the pass 0 butterfly */
73             
74             b = vec_perm(a1,a1,vcprmle(1,0,3,2));
75             b = vec_madd(a1,c1,b);
76             /* do the pass 0 butterfly */
77             
78             /* multiply third by -i */
79             b = vec_perm(b,b,vcprmle(2,3,1,0));
80             
81             /* do the pass 1 butterfly */
82             vec_st(vec_madd(b,c2,a), 0, r);
83             vec_st(vec_nmsub(b,c2,a), sizeof(vector float), r);
84             
85             r += 2;
86         } while (--j != 0);
87     }
88     /* pass 2 .. ln-1 */
89
90     nblocks = np >> 3;
91     nloops = 1 << 2;
92     np2 = np >> 1;
93
94     cptr1 = s->exptab1;
95     do {
96         p = z;
97         q = z + nloops;
98         j = nblocks;
99         do {
100             cptr = cptr1;
101             k = nloops >> 1;
102             do {
103                 vector float a,b,c,t1;
104
105                 a = vec_ld(0, (float*)p);
106                 b = vec_ld(0, (float*)q);
107                 
108                 /* complex mul */
109                 c = vec_ld(0, (float*)cptr);
110                 /*  cre*re cim*re */
111                 t1 = vec_madd(c, vec_perm(b,b,vcprmle(2,2,0,0)),vczero);
112                 c = vec_ld(sizeof(vector float), (float*)cptr);
113                 /*  -cim*im cre*im */
114                 b = vec_madd(c, vec_perm(b,b,vcprmle(3,3,1,1)),t1);
115                 
116                 /* butterfly */
117                 vec_st(vec_add(a,b), 0, (float*)p);
118                 vec_st(vec_sub(a,b), 0, (float*)q);
119                 
120                 p += 2;
121                 q += 2;
122                 cptr += 4;
123             } while (--k);
124             
125             p += nloops;
126             q += nloops;
127         } while (--j);
128         cptr1 += nloops * 2;
129         nblocks = nblocks >> 1;
130         nloops = nloops << 1;
131     } while (nblocks != 0);
132 }
133