altivec split-radix FFT (intrinsics version)
This is about 20% slower than the asm equivalent, due to function call and constant reloading overhead.
It's also from before I scheduled anything, doesn't have an unrolled fft16, and doesn't merge pass with swizzle.

diff --git a/libavcodec/fft.c b/libavcodec/fft.c
index a3f1151..d864f82 100644
--- a/libavcodec/fft.c
+++ b/libavcodec/fft.c
@@ -110,7 +110,6 @@ av_cold int ff_fft_init(FFTContext *s, int nbits, int inverse)
     has_vectors = mm_support();
     if (has_vectors & FF_MM_ALTIVEC) {
         s->fft_calc = ff_fft_calc_altivec;
-        split_radix = 0;
     }
 #endif
 
diff --git a/libavcodec/ppc/fft_altivec.c b/libavcodec/ppc/fft_altivec.c
index 7391131..d6463ae 100644
--- a/libavcodec/ppc/fft_altivec.c
+++ b/libavcodec/ppc/fft_altivec.c
@@ -1,8 +1,7 @@
 /*
  * FFT/IFFT transforms
  * AltiVec-enabled
- * Copyright (c) 2003 Romain Dolbeau <romain@dolbeau.org>
- * Based on code Copyright (c) 2002 Fabrice Bellard
+ * Copyright (c) 2009 Loren Merritt
  *
  * This file is part of FFmpeg.
  *
@@ -30,106 +29,157 @@
  * AltiVec-enabled
  * This code assumes that the 'z' pointer is 16 bytes-aligned
  * It also assumes all FFTComplex are 8 bytes-aligned pair of float
- * The code is exactly the same as the SSE version, except
- * that successive MUL + ADD/SUB have been merged into
- * fused multiply-add ('vec_madd' in altivec)
  */
-void ff_fft_calc_altivec(FFTContext *s, FFTComplex *z)
-{
-POWERPC_PERF_DECLARE(altivec_fft_num, s->nbits >= 6);
-    register const vector float vczero = (const vector float)vec_splat_u32(0.);
-
-    int ln = s->nbits;
-    int j, np, np2;
-    int nblocks, nloops;
-    register FFTComplex *p, *q;
-    FFTComplex *cptr, *cptr1;
-    int k;
-
-POWERPC_PERF_START_COUNT(altivec_fft_num, s->nbits >= 6);
-
-    np = 1 << ln;
-
-    {
-        vector float *r, a, b, a1, c1, c2;
-
-        r = (vector float *)&z[0];
-
-        c1 = vcii(p,p,n,n);
 
-        if (s->inverse) {
-            c2 = vcii(p,p,n,p);
-        } else {
-            c2 = vcii(p,p,p,n);
-        }
-
-        j = (np >> 2);
-        do {
-            a = vec_ld(0, r);
-            a1 = vec_ld(sizeof(vector float), r);
-
-            b = vec_perm(a,a,vcprmle(1,0,3,2));
-            a = vec_madd(a,c1,b);
-            /* do the pass 0 butterfly */
-
-            b = vec_perm(a1,a1,vcprmle(1,0,3,2));
-            b = vec_madd(a1,c1,b);
-            /* do the pass 0 butterfly */
+typedef vector float v4f;
+
+extern FFTSample ff_cos_16[];
+extern FFTSample ff_cos_32[];
+extern FFTSample ff_cos_64[];
+extern FFTSample ff_cos_128[];
+extern FFTSample ff_cos_256[];
+extern FFTSample ff_cos_512[];
+extern FFTSample ff_cos_1024[];
+extern FFTSample ff_cos_2048[];
+extern FFTSample ff_cos_4096[];
+extern FFTSample ff_cos_8192[];
+extern FFTSample ff_cos_16384[];
+extern FFTSample ff_cos_32768[];
+extern FFTSample ff_cos_65536[];
+
+// convert from simd order to C order
+static void deswizzle(v4f *z, int n)
+{
+    int i;
+    n >>= 1;
+    for(i=0; i<n; i+=2) {
+        v4f re = z[i];
+        v4f im = z[i+1];
+        z[i]   = vec_mergeh(re, im);
+        z[i+1] = vec_mergel(re, im);
+    }
+}
 
-            /* multiply third by -i */
-            b = vec_perm(b,b,vcprmle(2,3,1,0));
+#define BF(x,y,a,b) {\
+    v4f x1 = a - b;\
+    v4f y1 = a + b;\
+    x = x1;\
+    y = y1;\
+}
 
-            /* do the pass 1 butterfly */
-            vec_st(vec_madd(b,c2,a), 0, r);
-            vec_st(vec_nmsub(b,c2,a), sizeof(vector float), r);
+#define BUTTERFLIES(a0r,a0i,a1r,a1i,a2r,a2i,a3r,a3i) {\
+    BF(t3, t5, t5, t1);\
+    BF(a2r, a0r, a0r, t5);\
+    BF(a3i, a1i, a1i, t3);\
+    BF(t4, t6, t2, t6);\
+    BF(a3r, a1r, a1r, t4);\
+    BF(a2i, a0i, a0i, t6);\
+}
 
-            r += 2;
-        } while (--j != 0);
-    }
-    /* pass 2 .. ln-1 */
+#define TRANSFORM(a0r,a0i,a1r,a1i,a2r,a2i,a3r,a3i,wre,wim) {\
+    t1 = a2r * wre + a2i * wim;\
+    t2 = a2i * wre - a2r * wim;\
+    t5 = a3r * wre - a3i * wim;\
+    t6 = a3i * wre + a3r * wim;\
+    BUTTERFLIES(a0r,a0i,a1r,a1i,a2r,a2i,a3r,a3i)\
+}
 
-    nblocks = np >> 3;
-    nloops = 1 << 2;
-    np2 = np >> 1;
+/* z[0...8n-1], w[1...2n-1] */
+static void pass(v4f *z, const FFTSample *wre, unsigned int n)
+{
+    v4f t1, t2, t3, t4, t5, t6;
+    int o1 = 2*n;
+    int o2 = 4*n;
+    int o3 = 6*n;
+    const FFTSample *wim = wre+4*n;
 
-    cptr1 = s->exptab1;
     do {
-        p = z;
-        q = z + nloops;
-        j = nblocks;
-        do {
-            cptr = cptr1;
-            k = nloops >> 1;
-            do {
-                vector float a,b,c,t1;
-
-                a = vec_ld(0, (float*)p);
-                b = vec_ld(0, (float*)q);
+        v4f vwim = vec_perm(*(const v4f*)(wim-4), *(const v4f*)wim, vcprm(s0,3,2,1));
+        TRANSFORM(z[0], z[1], z[o1], z[o1+1], z[o2], z[o2+1], z[o3], z[o3+1], *(const v4f*)wre, vwim);
+        z += 2;
+        wre += 4;
+        wim -= 4;
+    } while(--n);
+}
 
-                /* complex mul */
-                c = vec_ld(0, (float*)cptr);
-                /*  cre*re cim*re */
-                t1 = vec_madd(c, vec_perm(b,b,vcprmle(2,2,0,0)),vczero);
-                c = vec_ld(sizeof(vector float), (float*)cptr);
-                /*  -cim*im cre*im */
-                b = vec_madd(c, vec_perm(b,b,vcprmle(3,3,1,1)),t1);
+#define DECL_FFT(n,n2,n4)\
+static void fft##n(v4f *z)\
+{\
+    fft##n2(z);\
+    fft##n4(z+n4);\
+    fft##n4(z+n4*3/2);\
+    pass(z,ff_cos_##n,n4/4);\
+}
 
-                /* butterfly */
-                vec_st(vec_add(a,b), 0, (float*)p);
-                vec_st(vec_sub(a,b), 0, (float*)q);
+// FIXME each different vec_perm costs 3 extra instructions: lis,addi,lvx.
+// some of the vec_perms can be changed to a vec_merge,
+// and the constants could be moved to a contiguous array to eliminate the address computation part.
+// better yet, load the constants into registers once and for all. but gcc isn't smart enough to keep them across function boundaries.
+static void fft4(v4f *z)
+{
+    v4f t1, t2, t3, t4;
+    t1 = vec_perm(z[0], z[1], vcprm(0,1,s2,s1)); // {r0,i0,r3,i2}
+    t2 = vec_perm(z[0], z[1], vcprm(2,3,s0,s3)); // {r1,i1,r2,i3}
+    t3 = t1 + t2; // {t1,t2,t6,t5}
+    t4 = t1 - t2; // {t3,t4,t8,t7}
+    t1 = vec_mergeh(t3, t4); // {t1,t3,t2,t4}
+    t2 = vec_perm(t3, t4, vcprm(2,s3,3,s2)); // {t6,t7,t5,t8}
+    t3 = t1 + t2; // {r0,r1,i0,i1}
+    t4 = t1 - t2; // {r2,r3,i2,i3}
+    z[0] = vec_perm(t3, t4, vcprm(0,1,s0,s1)); // {r0,r1,r2,r3}
+    z[1] = vec_perm(t3, t4, vcprm(2,3,s2,s3)); // {i0,i1,i2,i3}
+}
 
-                p += 2;
-                q += 2;
-                cptr += 4;
-            } while (--k);
+static void fft8(v4f *z)
+{
+    v4f t1, t2, t3, t4, t5, t6;
+    fft4(z);
+    t1 = vec_mergeh(z[2], z[3]); // {r4,r6,i4,i6}
+    t2 = vec_mergel(z[2], z[3]); // {r5,r7,i5,i7}
+    t3 = t1 + t2; // {t1,t3,t2,t4}
+    t4 = t1 - t2; // {r5,r7,i5,i7}
+    t5 = vec_perm(t4, t4, vcprm(2,3,0,1));
+    t4 *= (const v4f){-M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2}; // {-r5,r7,i5,-i7}
+    t5 *= (const v4f){ M_SQRT1_2, M_SQRT1_2, M_SQRT1_2,  M_SQRT1_2}; // {i5,i7,r5,r7}
+    t4 += t5; // {t8,ta,t7,t9}
+    t1 = vec_perm(t3, t4, vcprm(1,2,s3,s0)); // {t3,t2,t9,t8}
+    t2 = vec_perm(t3, t4, vcprm(0,3,s2,s1)); // {t1,t4,t7,ta}
+    t3 = t1 + t2; // {t1,t2,t9,ta}
+    t4 = t1 - t2; // {t6,t5,tc,tb}
+    t1 = vec_perm(t3, t4, vcprm(0,2,s1,s3)); // {t1,t9,t5,tb}
+    t2 = vec_perm(t3, t4, vcprm(1,3,s0,s2)); // {t2,ta,t6,tc}
+    t3 = z[0] + t1; // {r0,r1,r2,r3}
+    t4 = z[1] + t2; // {i0,i1,i2,i3}
+    t5 = z[0] - t1; // {r4,r5,r6,r7}
+    t6 = z[1] - t2; // {i4,i5,i6,i7}
+    z[0] = t3;
+    z[1] = t4;
+    z[2] = t5;
+    z[3] = t6;
+}
 
-            p += nloops;
-            q += nloops;
-        } while (--j);
-        cptr1 += nloops * 2;
-        nblocks = nblocks >> 1;
-        nloops = nloops << 1;
-    } while (nblocks != 0);
+DECL_FFT(16,8,4)
+DECL_FFT(32,16,8)
+DECL_FFT(64,32,16)
+DECL_FFT(128,64,32)
+DECL_FFT(256,128,64)
+DECL_FFT(512,256,128)
+DECL_FFT(1024,512,256)
+DECL_FFT(2048,1024,512)
+DECL_FFT(4096,2048,1024)
+DECL_FFT(8192,4096,2048)
+DECL_FFT(16384,8192,4096)
+DECL_FFT(32768,16384,8192)
+DECL_FFT(65536,32768,16384)
+
+static void (*fft_dispatch[])(v4f*) = {
+    fft4, fft8, fft16, fft32, fft64, fft128, fft256, fft512, fft1024,
+    fft2048, fft4096, fft8192, fft16384, fft32768, fft65536,
+};
 
-POWERPC_PERF_STOP_COUNT(altivec_fft_num, s->nbits >= 6);
+void ff_fft_calc_altivec(FFTContext *s, FFTComplex *z)
+{
+    fft_dispatch[s->nbits-2]((v4f*)z);
+    deswizzle((v4f*)z, 1<<s->nbits);
 }
+