[FFmpeg-devel] [PATCH] mdct15: add assembly optimizations for the 15-point FFT

Rostislav Pehlivanov atomnuker at gmail.com
Thu Jun 22 01:33:25 EEST 2017


fft15_c:    1961 decicycles in fft15, 4193738 runs,    566 skips
fft15_fma3:  995 decicycles in fft15, 4193984 runs,    320 skips

I believe the FFT5 macro can be improved further, so suggestions
are welcome.
---
 libavcodec/mdct15.c          | 186 +++++++++++++++++++++----------------------
 libavcodec/mdct15.h          |  26 +++---
 libavcodec/x86/Makefile      |   2 +
 libavcodec/x86/mdct15.asm    | 149 ++++++++++++++++++++++++++++++++++
 libavcodec/x86/mdct15_init.c |  93 ++++++++++++++++++++++
 5 files changed, 345 insertions(+), 111 deletions(-)
 create mode 100644 libavcodec/x86/mdct15.asm
 create mode 100644 libavcodec/x86/mdct15_init.c

diff --git a/libavcodec/mdct15.c b/libavcodec/mdct15.c
index 8c42ece483..0cf734caf9 100644
--- a/libavcodec/mdct15.c
+++ b/libavcodec/mdct15.c
@@ -57,11 +57,6 @@ av_cold void ff_mdct15_uninit(MDCT15Context **ps)
     av_freep(ps);
 }
 
-static void mdct15(MDCT15Context *s, float *dst, const float *src, ptrdiff_t stride);
-
-static void imdct15_half(MDCT15Context *s, float *dst, const float *src,
-                         ptrdiff_t stride, float scale);
-
 static inline int init_pfa_reindex_tabs(MDCT15Context *s)
 {
     int i, j;
@@ -93,88 +88,8 @@ static inline int init_pfa_reindex_tabs(MDCT15Context *s)
     return 0;
 }
 
-av_cold int ff_mdct15_init(MDCT15Context **ps, int inverse, int N, double scale)
-{
-    MDCT15Context *s;
-    double alpha, theta;
-    int len2 = 15 * (1 << N);
-    int len  = 2 * len2;
-    int i;
-
-    /* Tested and verified to work on everything in between */
-    if ((N < 2) || (N > 13))
-        return AVERROR(EINVAL);
-
-    s = av_mallocz(sizeof(*s));
-    if (!s)
-        return AVERROR(ENOMEM);
-
-    s->fft_n = N - 1;
-    s->len4 = len2 / 2;
-    s->len2 = len2;
-    s->inverse = inverse;
-    s->mdct = mdct15;
-    s->imdct_half = imdct15_half;
-
-    if (ff_fft_init(&s->ptwo_fft, N - 1, s->inverse) < 0)
-        goto fail;
-
-    if (init_pfa_reindex_tabs(s))
-        goto fail;
-
-    s->tmp  = av_malloc_array(len, 2 * sizeof(*s->tmp));
-    if (!s->tmp)
-        goto fail;
-
-    s->twiddle_exptab  = av_malloc_array(s->len4, sizeof(*s->twiddle_exptab));
-    if (!s->twiddle_exptab)
-        goto fail;
-
-    theta = 0.125f + (scale < 0 ? s->len4 : 0);
-    scale = sqrt(fabs(scale));
-    for (i = 0; i < s->len4; i++) {
-        alpha = 2 * M_PI * (i + theta) / len;
-        s->twiddle_exptab[i].re = cos(alpha) * scale;
-        s->twiddle_exptab[i].im = sin(alpha) * scale;
-    }
-
-    /* 15-point FFT exptab */
-    for (i = 0; i < 19; i++) {
-        if (i < 15) {
-            double theta = (2.0f * M_PI * i) / 15.0f;
-            if (!s->inverse)
-                theta *= -1;
-            s->exptab[i].re = cos(theta);
-            s->exptab[i].im = sin(theta);
-        } else { /* Wrap around to simplify fft15 */
-            s->exptab[i] = s->exptab[i - 15];
-        }
-    }
-
-    /* 5-point FFT exptab */
-    s->exptab[19].re = cos(2.0f * M_PI / 5.0f);
-    s->exptab[19].im = sin(2.0f * M_PI / 5.0f);
-    s->exptab[20].re = cos(1.0f * M_PI / 5.0f);
-    s->exptab[20].im = sin(1.0f * M_PI / 5.0f);
-
-    /* Invert the phase for an inverse transform, do nothing for a forward transform */
-    if (s->inverse) {
-        s->exptab[19].im *= -1;
-        s->exptab[20].im *= -1;
-    }
-
-    *ps = s;
-
-    return 0;
-
-fail:
-    ff_mdct15_uninit(&s);
-    return AVERROR(ENOMEM);
-}
-
 /* Stride is hardcoded to 3 */
-static inline void fft5(const FFTComplex exptab[2], FFTComplex *out,
-                        const FFTComplex *in)
+static inline void fft5(FFTComplex *out, FFTComplex *in, FFTComplex exptab[2])
 {
     FFTComplex z0[4], t[6];
 
@@ -219,14 +134,14 @@ static inline void fft5(const FFTComplex exptab[2], FFTComplex *out,
     out[4].im = in[0].im + z0[3].im;
 }
 
-static void fft15(const FFTComplex exptab[22], FFTComplex *out, const FFTComplex *in, size_t stride)
+static void fft15(FFTComplex *out, FFTComplex *in, FFTComplex *exptab, size_t stride)
 {
     int k;
     FFTComplex tmp1[5], tmp2[5], tmp3[5];
 
-    fft5(exptab + 19, tmp1, in + 0);
-    fft5(exptab + 19, tmp2, in + 1);
-    fft5(exptab + 19, tmp3, in + 2);
+    fft5(tmp1, in + 0, exptab + 19);
+    fft5(tmp2, in + 1, exptab + 19);
+    fft5(tmp3, in + 2, exptab + 19);
 
     for (k = 0; k < 5; k++) {
         FFTComplex t[2];
@@ -253,10 +168,10 @@ static void mdct15(MDCT15Context *s, float *dst, const float *src, ptrdiff_t str
     int i, j;
     const int len4 = s->len4, len3 = len4 * 3, len8 = len4 >> 1;
     const int l_ptwo = 1 << s->ptwo_fft.nbits;
-    FFTComplex fft15in[15];
 
     /* Folding and pre-reindexing */
     for (i = 0; i < l_ptwo; i++) {
+        DECLARE_ALIGNED(32, FFTComplex, fft15in)[32];
         for (j = 0; j < 15; j++) {
             float re, im;
             const int k = s->pfa_prereindex[i*15 + j];
@@ -269,7 +184,7 @@ static void mdct15(MDCT15Context *s, float *dst, const float *src, ptrdiff_t str
             }
             CMUL(fft15in[j].re, fft15in[j].im, re, im, s->twiddle_exptab[k].re, -s->twiddle_exptab[k].im);
         }
-        fft15(s->exptab, s->tmp + s->ptwo_fft.revtab[i], fft15in, l_ptwo);
+        s->fft15(s->tmp + s->ptwo_fft.revtab[i], fft15in, s->exptab, l_ptwo);
     }
 
     /* Then a 15xN FFT (where N is a power of two) */
@@ -294,19 +209,19 @@ static void mdct15(MDCT15Context *s, float *dst, const float *src, ptrdiff_t str
 static void imdct15_half(MDCT15Context *s, float *dst, const float *src,
                          ptrdiff_t stride, float scale)
 {
-    FFTComplex fft15in[15];
     FFTComplex *z = (FFTComplex *)dst;
     int i, j, len8 = s->len4 >> 1, l_ptwo = 1 << s->ptwo_fft.nbits;
     const float *in1 = src, *in2 = src + (s->len2 - 1) * stride;
 
     /* Reindex input, putting it into a buffer and doing an Nx15 FFT */
     for (i = 0; i < l_ptwo; i++) {
+        DECLARE_ALIGNED(32, FFTComplex, fft15in)[32];
         for (j = 0; j < 15; j++) {
             const int k = s->pfa_prereindex[i*15 + j];
             FFTComplex tmp = { *(in2 - 2*k*stride), *(in1 + 2*k*stride) };
             CMUL3(fft15in[j], tmp, s->twiddle_exptab[k]);
         }
-        fft15(s->exptab, s->tmp + s->ptwo_fft.revtab[i], fft15in, l_ptwo);
+        s->fft15(s->tmp + s->ptwo_fft.revtab[i], fft15in, s->exptab, l_ptwo);
     }
 
     /* Then a 15xN FFT (where N is a power of two) */
@@ -327,3 +242,86 @@ static void imdct15_half(MDCT15Context *s, float *dst, const float *src,
         z[i0].im = scale * im1;
     }
 }
+
+av_cold int ff_mdct15_init(MDCT15Context **ps, int inverse, int N, double scale)
+{
+    MDCT15Context *s;
+    double alpha, theta;
+    int len2 = 15 * (1 << N);
+    int len  = 2 * len2;
+    int i;
+
+    /* Tested and verified to work on everything in between */
+    if ((N < 2) || (N > 13))
+        return AVERROR(EINVAL);
+
+    s = av_mallocz(sizeof(*s));
+    if (!s)
+        return AVERROR(ENOMEM);
+
+    s->fft_n      = N - 1;
+    s->len4       = len2 / 2;
+    s->len2       = len2;
+    s->inverse    = inverse;
+    s->fft15      = fft15;
+    s->mdct       = mdct15;
+    s->imdct_half = imdct15_half;
+
+    if (ff_fft_init(&s->ptwo_fft, N - 1, s->inverse) < 0)
+        goto fail;
+
+    if (init_pfa_reindex_tabs(s))
+        goto fail;
+
+    s->tmp  = av_malloc_array(len, 2 * sizeof(*s->tmp));
+    if (!s->tmp)
+        goto fail;
+
+    s->twiddle_exptab  = av_malloc_array(s->len4, sizeof(*s->twiddle_exptab));
+    if (!s->twiddle_exptab)
+        goto fail;
+
+    theta = 0.125f + (scale < 0 ? s->len4 : 0);
+    scale = sqrt(fabs(scale));
+    for (i = 0; i < s->len4; i++) {
+        alpha = 2 * M_PI * (i + theta) / len;
+        s->twiddle_exptab[i].re = cos(alpha) * scale;
+        s->twiddle_exptab[i].im = sin(alpha) * scale;
+    }
+
+    /* 15-point FFT exptab */
+    for (i = 0; i < 19; i++) {
+        if (i < 15) {
+            double theta = (2.0f * M_PI * i) / 15.0f;
+            if (!s->inverse)
+                theta *= -1;
+            s->exptab[i].re = cos(theta);
+            s->exptab[i].im = sin(theta);
+        } else { /* Wrap around to simplify fft15 */
+            s->exptab[i] = s->exptab[i - 15];
+        }
+    }
+
+    /* 5-point FFT exptab */
+    s->exptab[19].re = cos(2.0f * M_PI / 5.0f);
+    s->exptab[19].im = sin(2.0f * M_PI / 5.0f);
+    s->exptab[20].re = cos(1.0f * M_PI / 5.0f);
+    s->exptab[20].im = sin(1.0f * M_PI / 5.0f);
+
+    /* Invert the phase for an inverse transform, do nothing for a forward transform */
+    if (s->inverse) {
+        s->exptab[19].im *= -1;
+        s->exptab[20].im *= -1;
+    }
+
+    if (ARCH_X86)
+        ff_mdct15_init_x86(s);
+
+    *ps = s;
+
+    return 0;
+
+fail:
+    ff_mdct15_uninit(&s);
+    return AVERROR(ENOMEM);
+}
diff --git a/libavcodec/mdct15.h b/libavcodec/mdct15.h
index ef94edff6c..15a91220b6 100644
--- a/libavcodec/mdct15.h
+++ b/libavcodec/mdct15.h
@@ -34,34 +34,26 @@ typedef struct MDCT15Context {
     int *pfa_postreindex;
 
     FFTContext ptwo_fft;
-
     FFTComplex *tmp;
-
     FFTComplex *twiddle_exptab;
 
-    /* 0 - 18: fft15 twiddles, 19 - 20: fft5 twiddles */
-    FFTComplex exptab[21];
+    DECLARE_ALIGNED(32, FFTComplex, exptab)[64];
 
-    /**
-     * Calculate a full 2N -> N MDCT
-     */
+    /* 15-point FFT */
+    void (*fft15)(FFTComplex *out, FFTComplex *in, FFTComplex *exptab, size_t stride);
+
+    /* Calculate a full 2N -> N MDCT */
     void (*mdct)(struct MDCT15Context *s, float *dst, const float *src, ptrdiff_t stride);
 
-    /**
-     * Calculate the middle half of the iMDCT
-     */
+    /* Calculate the middle half of the iMDCT */
     void (*imdct_half)(struct MDCT15Context *s, float *dst, const float *src,
                        ptrdiff_t src_stride, float scale);
 } MDCT15Context;
 
-/**
- * Init an (i)MDCT of the length 2 * 15 * (2^N)
- */
+/* Init an (i)MDCT of the length 2 * 15 * (2^N) */
 int ff_mdct15_init(MDCT15Context **ps, int inverse, int N, double scale);
-
-/**
- * Frees a context
- */
 void ff_mdct15_uninit(MDCT15Context **ps);
 
+void ff_mdct15_init_x86(MDCT15Context *s);
+
 #endif /* AVCODEC_MDCT15_H */
diff --git a/libavcodec/x86/Makefile b/libavcodec/x86/Makefile
index e65118d134..12f6b84372 100644
--- a/libavcodec/x86/Makefile
+++ b/libavcodec/x86/Makefile
@@ -53,6 +53,7 @@ OBJS-$(CONFIG_DCA_DECODER)             += x86/dcadsp_init.o x86/synth_filter_ini
 OBJS-$(CONFIG_DNXHD_ENCODER)           += x86/dnxhdenc_init.o
 OBJS-$(CONFIG_HEVC_DECODER)            += x86/hevcdsp_init.o
 OBJS-$(CONFIG_JPEG2000_DECODER)        += x86/jpeg2000dsp_init.o
+OBJS-$(CONFIG_MDCT15)                  += x86/mdct15_init.o
 OBJS-$(CONFIG_MLP_DECODER)             += x86/mlpdsp_init.o
 OBJS-$(CONFIG_MPEG4_DECODER)           += x86/xvididct_init.o
 OBJS-$(CONFIG_PNG_DECODER)             += x86/pngdsp_init.o
@@ -158,6 +159,7 @@ X86ASM-OBJS-$(CONFIG_HEVC_DECODER)     += x86/hevc_add_res.o            \
                                           x86/hevc_sao.o                \
                                           x86/hevc_sao_10bit.o
 X86ASM-OBJS-$(CONFIG_JPEG2000_DECODER) += x86/jpeg2000dsp.o
+X86ASM-OBJS-$(CONFIG_MDCT15)           += x86/mdct15.o
 X86ASM-OBJS-$(CONFIG_MLP_DECODER)      += x86/mlpdsp.o
 X86ASM-OBJS-$(CONFIG_MPEG4_DECODER)    += x86/xvididct.o
 X86ASM-OBJS-$(CONFIG_PNG_DECODER)      += x86/pngdsp.o
diff --git a/libavcodec/x86/mdct15.asm b/libavcodec/x86/mdct15.asm
new file mode 100644
index 0000000000..437df7bea6
--- /dev/null
+++ b/libavcodec/x86/mdct15.asm
@@ -0,0 +1,149 @@
+;******************************************************************************
+;* SIMD optimized non-power of two MDCT functions
+;*
+;* Copyright (C) 2017 Rostislav Pehlivanov <atomnuker at gmail.com>
+;*
+;* This file is part of FFmpeg.
+;*
+;* FFmpeg is free software; you can redistribute it and/or
+;* modify it under the terms of the GNU Lesser General Public
+;* License as published by the Free Software Foundation; either
+;* version 2.1 of the License, or (at your option) any later version.
+;*
+;* FFmpeg is distributed in the hope that it will be useful,
+;* but WITHOUT ANY WARRANTY; without even the implied warranty of
+;* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+;* Lesser General Public License for more details.
+;*
+;* You should have received a copy of the GNU Lesser General Public
+;* License along with FFmpeg; if not, write to the Free Software
+;* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+;******************************************************************************
+
+%include "libavutil/x86/x86util.asm"
+
+SECTION_RODATA 32
+
+sign_adjust_5: times 2 dd 1.0, -1.0, -1.0, 1.0
+
+SECTION .text
+
+%macro FFT5 7 ; src, offset, dst1 (1/2 filled), dst2, exptab1, exptab2, signadjust (uses m0-m4)
+    movddup xm0, [%1q + 0*16 + 0 + %2]  ; in[ 0].re, in[ 0].im, in[ 1].re, in[ 1].im
+    movsd   xm1, [%1q + 1*16 + 8 + %2]  ; in[ 3].re, in[ 3].im, in[ 4].re, in[ 4].im
+    movups  xm2, [%1q + 2*16 + 8 + %2]  ; in[ 5].re, in[ 5].im, in[ 6].re, in[ 6].im
+    movups  xm3, [%1q + 4*16 + 0 + %2]  ; in[ 8].re, in[ 8].im, in[ 9].re, in[ 9].im
+    movsd   xm4, [%1q + 6*16 + 0 + %2]  ; in[12].re, in[12].im, in[13].re, in[13].im
+
+    shufps      xm1,  xm2, q3210        ; in[ 3].re, in[ 3].im, in[ 6].re, in[ 6].im
+    shufps      xm4,  xm3, q3210        ; in[12].re, in[12].im, in[ 9].re, in[ 9].im
+
+    subps       xm2,  xm1, xm4          ; t[2].im, t[2].re, t[3].im, t[3].re
+    addps       xm1,  xm4               ; t[0].re, t[0].im, t[1].re, t[1].im
+
+    shufps      %3,   xm1, xm1, q3120   ; t[0].re, t[1].re, t[0].im, t[1].im
+    haddps      %3,   xm4
+    addps       %3,   xm0               ; DC[0].re, DC[0].im, junk...
+    movlhps     %3,   %3                ; DC[0].re, DC[0].im, DC[0].re, DC[0].im
+
+    shufps      xm3,  xm1, xm2, q0110   ; t[0].re, t[0].im, t[2].re, t[2].im
+    shufps      xm1,  xm2, q2332        ; t[1].re, t[1].im, t[3].re, t[3].im
+
+    mulps       xm%4, xm1, %5
+    mulps       xm4,  xm3, %6
+    mulps       xm1,  %6
+
+    mulps       xm1,  %7
+    fmaddsubps  xm3,  xm3, %5, xm1      ; t[0].re, t[0].im, t[2].re, t[2].im
+    subps       xm%4, xm4               ; t[4].re, t[4].im, t[5].re, t[5].im
+
+    mova        xm4,  %7
+
+    unpckhpd    xm2,  xm3, xm%4         ; t[2].re, t[2].im, t[5].re, t[5].im
+    unpcklpd    xm3,  xm%4              ; t[0].re, t[0].im, t[4].re, t[4].im
+
+    mova        xm%4, xm4
+    fmaddps     xm%4, xm2, xm%4, xm3
+    addps       xm%4, xm0               ; ACs(1)
+
+    fnmaddps    xm2,  xm2, xm4, xm3
+    shufps      xm2,  xm2, q1032
+    addps       xm2,  xm0               ; ACs(2)
+    vinsertf128 m%4,  m%4, xm2, 1       ; All ACs (tmp[1] through to tmp[4])
+%endmacro
+
+%macro BUTTERFLIES_DC 7 ; exptab, exptab_offset, src1, src2, src3, out, out_offset, (uses m0-m1)
+    movaps m0, [%1q + %2]
+    vextractf128 xm1, m0, 1
+
+    mulps xm1, %5
+    mulps xm0, %4
+
+    haddps xm0, xm1
+    shufps xm0, xm0, q3120              ; t[0].re, t[1].re, t[0].im, t[1].im
+
+    haddps xm0, xm0
+    addps  xm0, %3
+
+    movsd [%6q + %7], xm0
+%endmacro
+
+%macro BUTTERFLIES_AC 8 ; exptab, exptab_offset, src1, src2, src3, out, strideq, k_off
+    mulps m0, %4, [%1q + 64*0 + 0*mmsize + %2]
+    mulps m1, %4, [%1q + 64*0 + 1*mmsize + %2]
+    mulps m2, %5, [%1q + 64*1 + 0*mmsize + %2]
+    mulps m3, %5, [%1q + 64*1 + 1*mmsize + %2]
+
+    shufps m1, m1, q2301
+    shufps m3, m3, q2301
+
+    addps m0, m1
+    addps m2, m3
+    addps m0, m2
+    addps m0, %3
+
+    vextractf128 xmm1, m0, 1
+
+    mov r5q, %8
+    add r5q, strideq
+    movlps [%6q + r5q], xmm0
+
+    add r5q, strideq
+    movhps [%6q + r5q], xmm0
+
+    add r5q, strideq
+    movlps [%6q + r5q], xmm1
+
+    add r5q, strideq
+    movhps [%6q + r5q], xmm1
+%endmacro
+
+;***************************************************************************************
+;void ff_fft15_fma3(FFTComplex *out, FFTComplex *in, FFTComplex *exptab, size_t stride);
+;***************************************************************************************
+INIT_YMM fma3
+cglobal fft15, 4, 6, 14, out, in, exptab, stride
+
+    shl strideq, 3
+
+    movaps m5, [exptabq + 480]
+    vextractf128 xm6, m5, 1
+    movaps xm7, [sign_adjust_5]
+
+    FFT5 in,  0,  xm8, 11, xm5, xm6, xm7
+    FFT5 in,  8,  xm9, 12, xm5, xm6, xm7
+    FFT5 in, 16, xm10, 13, xm5, xm6, xm7
+
+    BUTTERFLIES_DC exptab, (8*6 + 4*0)*2*4, xm8, xm9, xm10, out,   0
+    imul r4q, strideq, 5
+    BUTTERFLIES_DC exptab, (8*6 + 4*1)*2*4, xm8, xm9, xm10, out, r4q
+    imul r4q, 2 ; stride * 10
+    BUTTERFLIES_DC exptab, (8*6 + 4*2)*2*4, xm8, xm9, xm10, out, r4q
+
+    BUTTERFLIES_AC exptab, (8*0)*2*4, m11, m12, m13, out, strideq,   0
+    imul r4q, strideq, 5
+    BUTTERFLIES_AC exptab, (8*2)*2*4, m11, m12, m13, out, strideq, r4q
+    imul r4q, 2 ; stride * 10
+    BUTTERFLIES_AC exptab, (8*4)*2*4, m11, m12, m13, out, strideq, r4q
+
+    RET
diff --git a/libavcodec/x86/mdct15_init.c b/libavcodec/x86/mdct15_init.c
new file mode 100644
index 0000000000..16306b893f
--- /dev/null
+++ b/libavcodec/x86/mdct15_init.c
@@ -0,0 +1,93 @@
+/*
+ * Non power of two MDCT assembly optimizations
+ * Copyright (C) 2017 Rostislav Pehlivanov <atomnuker at gmail.com>
+ *
+ * This file is part of FFmpeg.
+ *
+ * FFmpeg is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * FFmpeg is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with FFmpeg; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#include "config.h"
+
+#include "libavutil/float_dsp.h"
+#include "libavutil/x86/cpu.h"
+#include "libavcodec/mdct15.h"
+
+void ff_fft15_fma3(FFTComplex *out, FFTComplex *in, FFTComplex *exptab, size_t stride);
+
+static void perm_twiddles(MDCT15Context *s)
+{
+    int k;
+
+    FFTComplex tmp[64], tmp2[64];
+    memcpy(tmp, s->exptab, sizeof(FFTComplex)*32);
+
+    /* 15-point FFT twiddles */
+    for (k = 0; k < 5; k++) {
+        tmp2[6*k + 0] = tmp[k +  0];
+        tmp2[6*k + 2] = tmp[k +  5];
+        tmp2[6*k + 4] = tmp[k + 10];
+
+        tmp2[6*k + 1] = tmp[2 * (k + 0)];
+        tmp2[6*k + 3] = tmp[2 * (k + 5)];
+        tmp2[6*k + 5] = tmp[2 *  k + 5 ];
+    }
+
+    /* Newtab */
+
+    for (k = 0; k < 6; k++) {
+        FFTComplex ac_exp[] = {
+            { tmp2[6*1 + k].re,  tmp2[6*1 + k].re },
+            { tmp2[6*2 + k].re,  tmp2[6*2 + k].re },
+            { tmp2[6*3 + k].re,  tmp2[6*3 + k].re },
+            { tmp2[6*4 + k].re,  tmp2[6*4 + k].re },
+            { tmp2[6*1 + k].im, -tmp2[6*1 + k].im },
+            { tmp2[6*2 + k].im, -tmp2[6*2 + k].im },
+            { tmp2[6*3 + k].im, -tmp2[6*3 + k].im },
+            { tmp2[6*4 + k].im, -tmp2[6*4 + k].im },
+        };
+        memcpy(s->exptab + 8*k, ac_exp, 8*sizeof(FFTComplex));
+    }
+
+    /* Specialcase when k = 0 */
+    for (k = 0; k < 3; k++) {
+        FFTComplex dc_exp[] = {
+            { tmp2[2*k + 0].re, -tmp2[2*k + 0].im },
+            { tmp2[2*k + 0].im,  tmp2[2*k + 0].re },
+            { tmp2[2*k + 1].re, -tmp2[2*k + 1].im },
+            { tmp2[2*k + 1].im,  tmp2[2*k + 1].re },
+        };
+        memcpy(s->exptab + 8*6 + 4*k, dc_exp, 4*sizeof(FFTComplex));
+    }
+
+    /* 5-point FFT twiddles */
+    FFTComplex exp_5point[] = {
+        { tmp[19].re, tmp[19].re },
+        { tmp[19].im, tmp[19].im },
+        { tmp[20].re, tmp[20].re },
+        { tmp[20].im, tmp[20].im },
+    };
+    memcpy(s->exptab + 8*6 + 4*3, exp_5point, 4*sizeof(FFTComplex));
+}
+
+av_cold void ff_mdct15_init_x86(MDCT15Context *s)
+{
+    int cpu_flags = av_get_cpu_flags();
+
+    if (ARCH_X86_64 && EXTERNAL_FMA3(cpu_flags)) {
+        perm_twiddles(s);
+        s->fft15 = ff_fft15_fma3;
+    }
+}
-- 
2.13.1.611.g7e3b11ae1



More information about the ffmpeg-devel mailing list