FFmpeg coverage

GCC Code Coverage Report

Directory:	../../../ffmpeg/		Exec	Total	Coverage
File:	src/libavcodec/fft_template.c	Lines:	263	283	92.9 %
Date:	2020-09-25 23:16:12	Branches:	60	78	76.9 %


/*
 * FFT/IFFT transforms
 * Copyright (c) 2008 Loren Merritt
 * Copyright (c) 2002 Fabrice Bellard
 * Partly based on libdjbfft by D. J. Bernstein
 *
 * 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
 */

/**
 * @file
 * FFT/IFFT transforms.
 */

#include <stdlib.h>
#include <string.h>
#include "libavutil/mathematics.h"
#include "libavutil/thread.h"
#include "fft.h"
#include "fft-internal.h"

#if FFT_FIXED_32
#include "fft_table.h"

static void av_cold fft_lut_init(void)
{
    int n = 0;
    ff_fft_lut_init(ff_fft_offsets_lut, 0, 1 << 17, &n);
}

#else /* FFT_FIXED_32 */

/* cos(2*pi*x/n) for 0<=x<=n/4, followed by its reverse */
#if !CONFIG_HARDCODED_TABLES
COSTABLE(16);
COSTABLE(32);
COSTABLE(64);
COSTABLE(128);
COSTABLE(256);
COSTABLE(512);
COSTABLE(1024);
COSTABLE(2048);
COSTABLE(4096);
COSTABLE(8192);
COSTABLE(16384);
COSTABLE(32768);
COSTABLE(65536);
COSTABLE(131072);

static av_cold void init_ff_cos_tabs(int index)
{
    int i;
    int m = 1<<index;
    double freq = 2*M_PI/m;
    FFTSample *tab = FFT_NAME(ff_cos_tabs)[index];
    for(i=0; i<=m/4; i++)
        tab[i] = FIX15(cos(i*freq));
    for(i=1; i<m/4; i++)
        tab[m/2-i] = tab[i];
}

typedef struct CosTabsInitOnce {
    void (*func)(void);
    AVOnce control;
} CosTabsInitOnce;

#define INIT_FF_COS_TABS_FUNC(index, size)          \
static av_cold void init_ff_cos_tabs_ ## size (void)\
{                                                   \
    init_ff_cos_tabs(index);                        \
}

INIT_FF_COS_TABS_FUNC(4, 16)
INIT_FF_COS_TABS_FUNC(5, 32)
INIT_FF_COS_TABS_FUNC(6, 64)
INIT_FF_COS_TABS_FUNC(7, 128)
INIT_FF_COS_TABS_FUNC(8, 256)
INIT_FF_COS_TABS_FUNC(9, 512)
INIT_FF_COS_TABS_FUNC(10, 1024)
INIT_FF_COS_TABS_FUNC(11, 2048)
INIT_FF_COS_TABS_FUNC(12, 4096)
INIT_FF_COS_TABS_FUNC(13, 8192)
INIT_FF_COS_TABS_FUNC(14, 16384)
INIT_FF_COS_TABS_FUNC(15, 32768)
INIT_FF_COS_TABS_FUNC(16, 65536)
INIT_FF_COS_TABS_FUNC(17, 131072)

static CosTabsInitOnce cos_tabs_init_once[] = {
    { NULL },
    { NULL },
    { NULL },
    { NULL },
    { init_ff_cos_tabs_16, AV_ONCE_INIT },
    { init_ff_cos_tabs_32, AV_ONCE_INIT },
    { init_ff_cos_tabs_64, AV_ONCE_INIT },
    { init_ff_cos_tabs_128, AV_ONCE_INIT },
    { init_ff_cos_tabs_256, AV_ONCE_INIT },
    { init_ff_cos_tabs_512, AV_ONCE_INIT },
    { init_ff_cos_tabs_1024, AV_ONCE_INIT },
    { init_ff_cos_tabs_2048, AV_ONCE_INIT },
    { init_ff_cos_tabs_4096, AV_ONCE_INIT },
    { init_ff_cos_tabs_8192, AV_ONCE_INIT },
    { init_ff_cos_tabs_16384, AV_ONCE_INIT },
    { init_ff_cos_tabs_32768, AV_ONCE_INIT },
    { init_ff_cos_tabs_65536, AV_ONCE_INIT },
    { init_ff_cos_tabs_131072, AV_ONCE_INIT },
};

#endif
COSTABLE_CONST FFTSample * const FFT_NAME(ff_cos_tabs)[] = {
    NULL, NULL, NULL, NULL,
    FFT_NAME(ff_cos_16),
    FFT_NAME(ff_cos_32),
    FFT_NAME(ff_cos_64),
    FFT_NAME(ff_cos_128),
    FFT_NAME(ff_cos_256),
    FFT_NAME(ff_cos_512),
    FFT_NAME(ff_cos_1024),
    FFT_NAME(ff_cos_2048),
    FFT_NAME(ff_cos_4096),
    FFT_NAME(ff_cos_8192),
    FFT_NAME(ff_cos_16384),
    FFT_NAME(ff_cos_32768),
    FFT_NAME(ff_cos_65536),
    FFT_NAME(ff_cos_131072),
};

#endif /* FFT_FIXED_32 */

static void fft_permute_c(FFTContext *s, FFTComplex *z);
static void fft_calc_c(FFTContext *s, FFTComplex *z);

static int split_radix_permutation(int i, int n, int inverse)
{
    int m;
    if(n <= 2) return i&1;
    m = n >> 1;
    if(!(i&m))            return split_radix_permutation(i, m, inverse)*2;
    m >>= 1;
    if(inverse == !(i&m)) return split_radix_permutation(i, m, inverse)*4 + 1;
    else                  return split_radix_permutation(i, m, inverse)*4 - 1;
}

av_cold void ff_init_ff_cos_tabs(int index)
{
#if (!CONFIG_HARDCODED_TABLES) && (!FFT_FIXED_32)
    ff_thread_once(&cos_tabs_init_once[index].control, cos_tabs_init_once[index].func);
#endif
}

static const int avx_tab[] = {
    0, 4, 1, 5, 8, 12, 9, 13, 2, 6, 3, 7, 10, 14, 11, 15
};

static int is_second_half_of_fft32(int i, int n)
{
    if (n <= 32)
        return i >= 16;
    else if (i < n/2)
        return is_second_half_of_fft32(i, n/2);
    else if (i < 3*n/4)
        return is_second_half_of_fft32(i - n/2, n/4);
    else
        return is_second_half_of_fft32(i - 3*n/4, n/4);
}

static av_cold void fft_perm_avx(FFTContext *s)
{
    int i;
    int n = 1 << s->nbits;

    for (i = 0; i < n; i += 16) {
        int k;
        if (is_second_half_of_fft32(i, n)) {
            for (k = 0; k < 16; k++)
                s->revtab[-split_radix_permutation(i + k, n, s->inverse) & (n - 1)] =
                    i + avx_tab[k];

        } else {
            for (k = 0; k < 16; k++) {
                int j = i + k;
                j = (j & ~7) | ((j >> 1) & 3) | ((j << 2) & 4);
                s->revtab[-split_radix_permutation(i + k, n, s->inverse) & (n - 1)] = j;
            }
        }
    }
}

av_cold int ff_fft_init(FFTContext *s, int nbits, int inverse)
{
    int i, j, n;

    s->revtab = NULL;
    s->revtab32 = NULL;

    if (nbits < 2 || nbits > 17)
        goto fail;
    s->nbits = nbits;
    n = 1 << nbits;

    if (nbits <= 16) {
        s->revtab = av_malloc(n * sizeof(uint16_t));
        if (!s->revtab)
            goto fail;
    } else {
        s->revtab32 = av_malloc(n * sizeof(uint32_t));
        if (!s->revtab32)
            goto fail;
    }
    s->tmp_buf = av_malloc(n * sizeof(FFTComplex));
    if (!s->tmp_buf)
        goto fail;
    s->inverse = inverse;
    s->fft_permutation = FF_FFT_PERM_DEFAULT;

    s->fft_permute = fft_permute_c;
    s->fft_calc    = fft_calc_c;
#if CONFIG_MDCT
    s->imdct_calc  = ff_imdct_calc_c;
    s->imdct_half  = ff_imdct_half_c;
    s->mdct_calc   = ff_mdct_calc_c;
#endif

#if FFT_FIXED_32
    {
        static AVOnce control = AV_ONCE_INIT;
        ff_thread_once(&control, fft_lut_init);
    }
#else /* FFT_FIXED_32 */
#if FFT_FLOAT
    if (ARCH_AARCH64) ff_fft_init_aarch64(s);
    if (ARCH_ARM)     ff_fft_init_arm(s);
    if (ARCH_PPC)     ff_fft_init_ppc(s);
    if (ARCH_X86)     ff_fft_init_x86(s);
    if (CONFIG_MDCT)  s->mdct_calcw = s->mdct_calc;
    if (HAVE_MIPSFPU) ff_fft_init_mips(s);
#else
    if (CONFIG_MDCT)  s->mdct_calcw = ff_mdct_calcw_c;
    if (ARCH_ARM)     ff_fft_fixed_init_arm(s);
#endif
    for(j=4; j<=nbits; j++) {
        ff_init_ff_cos_tabs(j);
    }
#endif /* FFT_FIXED_32 */


    if (s->fft_permutation == FF_FFT_PERM_AVX) {
        fft_perm_avx(s);
    } else {
#define PROCESS_FFT_PERM_SWAP_LSBS(num) do {\
    for(i = 0; i < n; i++) {\
        int k;\
        j = i;\
        j = (j & ~3) | ((j >> 1) & 1) | ((j << 1) & 2);\
        k = -split_radix_permutation(i, n, s->inverse) & (n - 1);\
        s->revtab##num[k] = j;\
    } \
} while(0);

#define PROCESS_FFT_PERM_DEFAULT(num) do {\
    for(i = 0; i < n; i++) {\
        int k;\
        j = i;\
        k = -split_radix_permutation(i, n, s->inverse) & (n - 1);\
        s->revtab##num[k] = j;\
    } \
} while(0);

#define SPLIT_RADIX_PERMUTATION(num) do { \
    if (s->fft_permutation == FF_FFT_PERM_SWAP_LSBS) {\
        PROCESS_FFT_PERM_SWAP_LSBS(num) \
    } else {\
        PROCESS_FFT_PERM_DEFAULT(num) \
    }\
} while(0);

    if (s->revtab)
        SPLIT_RADIX_PERMUTATION()
    if (s->revtab32)
        SPLIT_RADIX_PERMUTATION(32)

#undef PROCESS_FFT_PERM_DEFAULT
#undef PROCESS_FFT_PERM_SWAP_LSBS
#undef SPLIT_RADIX_PERMUTATION
    }

    return 0;
 fail:
    av_freep(&s->revtab);
    av_freep(&s->revtab32);
    av_freep(&s->tmp_buf);
    return -1;
}

static void fft_permute_c(FFTContext *s, FFTComplex *z)
{
    int j, np;
    const uint16_t *revtab = s->revtab;
    const uint32_t *revtab32 = s->revtab32;
    np = 1 << s->nbits;
    /* TODO: handle split-radix permute in a more optimal way, probably in-place */
    if (revtab) {
        for(j=0;j<np;j++) s->tmp_buf[revtab[j]] = z[j];
    } else
        for(j=0;j<np;j++) s->tmp_buf[revtab32[j]] = z[j];

    memcpy(z, s->tmp_buf, np * sizeof(FFTComplex));
}

av_cold void ff_fft_end(FFTContext *s)
{
    av_freep(&s->revtab);
    av_freep(&s->revtab32);
    av_freep(&s->tmp_buf);
}

#if FFT_FIXED_32

static void fft_calc_c(FFTContext *s, FFTComplex *z) {

    int nbits, i, n, num_transforms, offset, step;
    int n4, n2, n34;
    unsigned tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8;
    FFTComplex *tmpz;
    const int fft_size = (1 << s->nbits);
    int64_t accu;

    num_transforms = (0x2aab >> (16 - s->nbits)) | 1;

    for (n=0; n<num_transforms; n++){
        offset = ff_fft_offsets_lut[n] << 2;
        tmpz = z + offset;

        tmp1 = tmpz[0].re + (unsigned)tmpz[1].re;
        tmp5 = tmpz[2].re + (unsigned)tmpz[3].re;
        tmp2 = tmpz[0].im + (unsigned)tmpz[1].im;
        tmp6 = tmpz[2].im + (unsigned)tmpz[3].im;
        tmp3 = tmpz[0].re - (unsigned)tmpz[1].re;
        tmp8 = tmpz[2].im - (unsigned)tmpz[3].im;
        tmp4 = tmpz[0].im - (unsigned)tmpz[1].im;
        tmp7 = tmpz[2].re - (unsigned)tmpz[3].re;

        tmpz[0].re = tmp1 + tmp5;
        tmpz[2].re = tmp1 - tmp5;
        tmpz[0].im = tmp2 + tmp6;
        tmpz[2].im = tmp2 - tmp6;
        tmpz[1].re = tmp3 + tmp8;
        tmpz[3].re = tmp3 - tmp8;
        tmpz[1].im = tmp4 - tmp7;
        tmpz[3].im = tmp4 + tmp7;
    }

    if (fft_size < 8)
        return;

    num_transforms = (num_transforms >> 1) | 1;

    for (n=0; n<num_transforms; n++){
        offset = ff_fft_offsets_lut[n] << 3;
        tmpz = z + offset;

        tmp1 = tmpz[4].re + (unsigned)tmpz[5].re;
        tmp3 = tmpz[6].re + (unsigned)tmpz[7].re;
        tmp2 = tmpz[4].im + (unsigned)tmpz[5].im;
        tmp4 = tmpz[6].im + (unsigned)tmpz[7].im;
        tmp5 = tmp1 + tmp3;
        tmp7 = tmp1 - tmp3;
        tmp6 = tmp2 + tmp4;
        tmp8 = tmp2 - tmp4;

        tmp1 = tmpz[4].re - (unsigned)tmpz[5].re;
        tmp2 = tmpz[4].im - (unsigned)tmpz[5].im;
        tmp3 = tmpz[6].re - (unsigned)tmpz[7].re;
        tmp4 = tmpz[6].im - (unsigned)tmpz[7].im;

        tmpz[4].re = tmpz[0].re - tmp5;
        tmpz[0].re = tmpz[0].re + tmp5;
        tmpz[4].im = tmpz[0].im - tmp6;
        tmpz[0].im = tmpz[0].im + tmp6;
        tmpz[6].re = tmpz[2].re - tmp8;
        tmpz[2].re = tmpz[2].re + tmp8;
        tmpz[6].im = tmpz[2].im + tmp7;
        tmpz[2].im = tmpz[2].im - tmp7;

        accu = (int64_t)Q31(M_SQRT1_2)*(int)(tmp1 + tmp2);
        tmp5 = (int32_t)((accu + 0x40000000) >> 31);
        accu = (int64_t)Q31(M_SQRT1_2)*(int)(tmp3 - tmp4);
        tmp7 = (int32_t)((accu + 0x40000000) >> 31);
        accu = (int64_t)Q31(M_SQRT1_2)*(int)(tmp2 - tmp1);
        tmp6 = (int32_t)((accu + 0x40000000) >> 31);
        accu = (int64_t)Q31(M_SQRT1_2)*(int)(tmp3 + tmp4);
        tmp8 = (int32_t)((accu + 0x40000000) >> 31);
        tmp1 = tmp5 + tmp7;
        tmp3 = tmp5 - tmp7;
        tmp2 = tmp6 + tmp8;
        tmp4 = tmp6 - tmp8;

        tmpz[5].re = tmpz[1].re - tmp1;
        tmpz[1].re = tmpz[1].re + tmp1;
        tmpz[5].im = tmpz[1].im - tmp2;
        tmpz[1].im = tmpz[1].im + tmp2;
        tmpz[7].re = tmpz[3].re - tmp4;
        tmpz[3].re = tmpz[3].re + tmp4;
        tmpz[7].im = tmpz[3].im + tmp3;
        tmpz[3].im = tmpz[3].im - tmp3;
    }

    step = 1 << ((MAX_LOG2_NFFT-4) - 4);
    n4 = 4;

    for (nbits=4; nbits<=s->nbits; nbits++){
        n2  = 2*n4;
        n34 = 3*n4;
        num_transforms = (num_transforms >> 1) | 1;

        for (n=0; n<num_transforms; n++){
            const FFTSample *w_re_ptr = ff_w_tab_sr + step;
            const FFTSample *w_im_ptr = ff_w_tab_sr + MAX_FFT_SIZE/(4*16) - step;
            offset = ff_fft_offsets_lut[n] << nbits;
            tmpz = z + offset;

            tmp5 = tmpz[ n2].re + (unsigned)tmpz[n34].re;
            tmp1 = tmpz[ n2].re - (unsigned)tmpz[n34].re;
            tmp6 = tmpz[ n2].im + (unsigned)tmpz[n34].im;
            tmp2 = tmpz[ n2].im - (unsigned)tmpz[n34].im;

            tmpz[ n2].re = tmpz[ 0].re - tmp5;
            tmpz[  0].re = tmpz[ 0].re + tmp5;
            tmpz[ n2].im = tmpz[ 0].im - tmp6;
            tmpz[  0].im = tmpz[ 0].im + tmp6;
            tmpz[n34].re = tmpz[n4].re - tmp2;
            tmpz[ n4].re = tmpz[n4].re + tmp2;
            tmpz[n34].im = tmpz[n4].im + tmp1;
            tmpz[ n4].im = tmpz[n4].im - tmp1;

            for (i=1; i<n4; i++){
                FFTSample w_re = w_re_ptr[0];
                FFTSample w_im = w_im_ptr[0];
                accu  = (int64_t)w_re*tmpz[ n2+i].re;
                accu += (int64_t)w_im*tmpz[ n2+i].im;
                tmp1 = (int32_t)((accu + 0x40000000) >> 31);
                accu  = (int64_t)w_re*tmpz[ n2+i].im;
                accu -= (int64_t)w_im*tmpz[ n2+i].re;
                tmp2 = (int32_t)((accu + 0x40000000) >> 31);
                accu  = (int64_t)w_re*tmpz[n34+i].re;
                accu -= (int64_t)w_im*tmpz[n34+i].im;
                tmp3 = (int32_t)((accu + 0x40000000) >> 31);
                accu  = (int64_t)w_re*tmpz[n34+i].im;
                accu += (int64_t)w_im*tmpz[n34+i].re;
                tmp4 = (int32_t)((accu + 0x40000000) >> 31);

                tmp5 = tmp1 + tmp3;
                tmp1 = tmp1 - tmp3;
                tmp6 = tmp2 + tmp4;
                tmp2 = tmp2 - tmp4;

                tmpz[ n2+i].re = tmpz[   i].re - tmp5;
                tmpz[    i].re = tmpz[   i].re + tmp5;
                tmpz[ n2+i].im = tmpz[   i].im - tmp6;
                tmpz[    i].im = tmpz[   i].im + tmp6;
                tmpz[n34+i].re = tmpz[n4+i].re - tmp2;
                tmpz[ n4+i].re = tmpz[n4+i].re + tmp2;
                tmpz[n34+i].im = tmpz[n4+i].im + tmp1;
                tmpz[ n4+i].im = tmpz[n4+i].im - tmp1;

                w_re_ptr += step;
                w_im_ptr -= step;
            }
        }
        step >>= 1;
        n4   <<= 1;
    }
}

#else /* FFT_FIXED_32 */

#define BUTTERFLIES(a0,a1,a2,a3) {\
    BF(t3, t5, t5, t1);\
    BF(a2.re, a0.re, a0.re, t5);\
    BF(a3.im, a1.im, a1.im, t3);\
    BF(t4, t6, t2, t6);\
    BF(a3.re, a1.re, a1.re, t4);\
    BF(a2.im, a0.im, a0.im, t6);\
}

// force loading all the inputs before storing any.
// this is slightly slower for small data, but avoids store->load aliasing
// for addresses separated by large powers of 2.
#define BUTTERFLIES_BIG(a0,a1,a2,a3) {\
    FFTSample r0=a0.re, i0=a0.im, r1=a1.re, i1=a1.im;\
    BF(t3, t5, t5, t1);\
    BF(a2.re, a0.re, r0, t5);\
    BF(a3.im, a1.im, i1, t3);\
    BF(t4, t6, t2, t6);\
    BF(a3.re, a1.re, r1, t4);\
    BF(a2.im, a0.im, i0, t6);\
}

#define TRANSFORM(a0,a1,a2,a3,wre,wim) {\
    CMUL(t1, t2, a2.re, a2.im, wre, -wim);\
    CMUL(t5, t6, a3.re, a3.im, wre,  wim);\
    BUTTERFLIES(a0,a1,a2,a3)\
}

#define TRANSFORM_ZERO(a0,a1,a2,a3) {\
    t1 = a2.re;\
    t2 = a2.im;\
    t5 = a3.re;\
    t6 = a3.im;\
    BUTTERFLIES(a0,a1,a2,a3)\
}

/* z[0...8n-1], w[1...2n-1] */
#define PASS(name)\
static void name(FFTComplex *z, const FFTSample *wre, unsigned int n)\
{\
    FFTDouble t1, t2, t3, t4, t5, t6;\
    int o1 = 2*n;\
    int o2 = 4*n;\
    int o3 = 6*n;\
    const FFTSample *wim = wre+o1;\
    n--;\
\
    TRANSFORM_ZERO(z[0],z[o1],z[o2],z[o3]);\
    TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
    do {\
        z += 2;\
        wre += 2;\
        wim -= 2;\
        TRANSFORM(z[0],z[o1],z[o2],z[o3],wre[0],wim[0]);\
        TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
    } while(--n);\
}

PASS(pass)
#if !CONFIG_SMALL
#undef BUTTERFLIES
#define BUTTERFLIES BUTTERFLIES_BIG
PASS(pass_big)
#endif

#define DECL_FFT(n,n2,n4)\
static void fft##n(FFTComplex *z)\
{\
    fft##n2(z);\
    fft##n4(z+n4*2);\
    fft##n4(z+n4*3);\
    pass(z,FFT_NAME(ff_cos_##n),n4/2);\
}

static void fft4(FFTComplex *z)
{
    FFTDouble t1, t2, t3, t4, t5, t6, t7, t8;

    BF(t3, t1, z[0].re, z[1].re);
    BF(t8, t6, z[3].re, z[2].re);
    BF(z[2].re, z[0].re, t1, t6);
    BF(t4, t2, z[0].im, z[1].im);
    BF(t7, t5, z[2].im, z[3].im);
    BF(z[3].im, z[1].im, t4, t8);
    BF(z[3].re, z[1].re, t3, t7);
    BF(z[2].im, z[0].im, t2, t5);
}

static void fft8(FFTComplex *z)
{
    FFTDouble t1, t2, t3, t4, t5, t6;

    fft4(z);

    BF(t1, z[5].re, z[4].re, -z[5].re);
    BF(t2, z[5].im, z[4].im, -z[5].im);
    BF(t5, z[7].re, z[6].re, -z[7].re);
    BF(t6, z[7].im, z[6].im, -z[7].im);

    BUTTERFLIES(z[0],z[2],z[4],z[6]);
    TRANSFORM(z[1],z[3],z[5],z[7],sqrthalf,sqrthalf);
}

#if !CONFIG_SMALL
static void fft16(FFTComplex *z)
{
    FFTDouble t1, t2, t3, t4, t5, t6;
    FFTSample cos_16_1 = FFT_NAME(ff_cos_16)[1];
    FFTSample cos_16_3 = FFT_NAME(ff_cos_16)[3];

    fft8(z);
    fft4(z+8);
    fft4(z+12);

    TRANSFORM_ZERO(z[0],z[4],z[8],z[12]);
    TRANSFORM(z[2],z[6],z[10],z[14],sqrthalf,sqrthalf);
    TRANSFORM(z[1],z[5],z[9],z[13],cos_16_1,cos_16_3);
    TRANSFORM(z[3],z[7],z[11],z[15],cos_16_3,cos_16_1);
}
#else
DECL_FFT(16,8,4)
#endif
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)
#if !CONFIG_SMALL
#define pass pass_big
#endif
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)
DECL_FFT(131072,65536,32768)

static void (* const fft_dispatch[])(FFTComplex*) = {
    fft4, fft8, fft16, fft32, fft64, fft128, fft256, fft512, fft1024,
    fft2048, fft4096, fft8192, fft16384, fft32768, fft65536, fft131072
};

static void fft_calc_c(FFTContext *s, FFTComplex *z)
{
    fft_dispatch[s->nbits-2](z);
}
#endif /* FFT_FIXED_32 */


Generated by: GCOVR (Version 4.2)

Line	Branch	Exec	Source
1			/*
2			* FFT/IFFT transforms
3			* Copyright (c) 2008 Loren Merritt
4			* Copyright (c) 2002 Fabrice Bellard
5			* Partly based on libdjbfft by D. J. Bernstein
6			*
7			* This file is part of FFmpeg.
8			*
9			* FFmpeg is free software; you can redistribute it and/or
10			* modify it under the terms of the GNU Lesser General Public
11			* License as published by the Free Software Foundation; either
12			* version 2.1 of the License, or (at your option) any later version.
13			*
14			* FFmpeg is distributed in the hope that it will be useful,
15			* but WITHOUT ANY WARRANTY; without even the implied warranty of
16			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17			* Lesser General Public License for more details.
18			*
19			* You should have received a copy of the GNU Lesser General Public
20			* License along with FFmpeg; if not, write to the Free Software
21			* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22			*/
23
24			/**
25			* @file
26			* FFT/IFFT transforms.
27			*/
28
29			#include <stdlib.h>
30			#include <string.h>
31			#include "libavutil/mathematics.h"
32			#include "libavutil/thread.h"
33			#include "fft.h"
34			#include "fft-internal.h"
35
36			#if FFT_FIXED_32
37			#include "fft_table.h"
38
39		58	static void av_cold fft_lut_init(void)
40			{
41		58	int n = 0;
42		58	ff_fft_lut_init(ff_fft_offsets_lut, 0, 1 << 17, &n);
43		58	}
44
45			#else /* FFT_FIXED_32 */
46
47			/* cos(2pix/n) for 0<=x<=n/4, followed by its reverse */
48			#if !CONFIG_HARDCODED_TABLES
49			COSTABLE(16);
50			COSTABLE(32);
51			COSTABLE(64);
52			COSTABLE(128);
53			COSTABLE(256);
54			COSTABLE(512);
55			COSTABLE(1024);
56			COSTABLE(2048);
57			COSTABLE(4096);
58			COSTABLE(8192);
59			COSTABLE(16384);
60			COSTABLE(32768);
61			COSTABLE(65536);
62			COSTABLE(131072);
63
64		2588	static av_cold void init_ff_cos_tabs(int index)
65			{
66			int i;
67		2588	int m = 1<<index;
68		2588	double freq = 2*M_PI/m;
69		2588	FFTSample *tab = FFT_NAME(ff_cos_tabs)[index];
70	✓✓	191096	for(i=0; i<=m/4; i++)
71		188508	tab[i] = FIX15(cos(i*freq));
72	✓✓	185920	for(i=1; i<m/4; i++)
73		183332	tab[m/2-i] = tab[i];
74		2588	}
75
76			typedef struct CosTabsInitOnce {
77			void (*func)(void);
78			AVOnce control;
79			} CosTabsInitOnce;
80
81			#define INIT_FF_COS_TABS_FUNC(index, size) \
82			static av_cold void init_ff_cos_tabs_ ## size (void)\
83			{ \
84			init_ff_cos_tabs(index); \
85			}
86
87		514	INIT_FF_COS_TABS_FUNC(4, 16)
88		493	INIT_FF_COS_TABS_FUNC(5, 32)
89		415	INIT_FF_COS_TABS_FUNC(6, 64)
90		387	INIT_FF_COS_TABS_FUNC(7, 128)
91		326	INIT_FF_COS_TABS_FUNC(8, 256)
92		298	INIT_FF_COS_TABS_FUNC(9, 512)
93		76	INIT_FF_COS_TABS_FUNC(10, 1024)
94		41	INIT_FF_COS_TABS_FUNC(11, 2048)
95		26	INIT_FF_COS_TABS_FUNC(12, 4096)
96		7	INIT_FF_COS_TABS_FUNC(13, 8192)
97		5	INIT_FF_COS_TABS_FUNC(14, 16384)
98			INIT_FF_COS_TABS_FUNC(15, 32768)
99			INIT_FF_COS_TABS_FUNC(16, 65536)
100			INIT_FF_COS_TABS_FUNC(17, 131072)
101
102			static CosTabsInitOnce cos_tabs_init_once[] = {
103			{ NULL },
104			{ NULL },
105			{ NULL },
106			{ NULL },
107			{ init_ff_cos_tabs_16, AV_ONCE_INIT },
108			{ init_ff_cos_tabs_32, AV_ONCE_INIT },
109			{ init_ff_cos_tabs_64, AV_ONCE_INIT },
110			{ init_ff_cos_tabs_128, AV_ONCE_INIT },
111			{ init_ff_cos_tabs_256, AV_ONCE_INIT },
112			{ init_ff_cos_tabs_512, AV_ONCE_INIT },
113			{ init_ff_cos_tabs_1024, AV_ONCE_INIT },
114			{ init_ff_cos_tabs_2048, AV_ONCE_INIT },
115			{ init_ff_cos_tabs_4096, AV_ONCE_INIT },
116			{ init_ff_cos_tabs_8192, AV_ONCE_INIT },
117			{ init_ff_cos_tabs_16384, AV_ONCE_INIT },
118			{ init_ff_cos_tabs_32768, AV_ONCE_INIT },
119			{ init_ff_cos_tabs_65536, AV_ONCE_INIT },
120			{ init_ff_cos_tabs_131072, AV_ONCE_INIT },
121			};
122
123			#endif
124			COSTABLE_CONST FFTSample * const FFT_NAME(ff_cos_tabs)[] = {
125			NULL, NULL, NULL, NULL,
126			FFT_NAME(ff_cos_16),
127			FFT_NAME(ff_cos_32),
128			FFT_NAME(ff_cos_64),
129			FFT_NAME(ff_cos_128),
130			FFT_NAME(ff_cos_256),
131			FFT_NAME(ff_cos_512),
132			FFT_NAME(ff_cos_1024),
133			FFT_NAME(ff_cos_2048),
134			FFT_NAME(ff_cos_4096),
135			FFT_NAME(ff_cos_8192),
136			FFT_NAME(ff_cos_16384),
137			FFT_NAME(ff_cos_32768),
138			FFT_NAME(ff_cos_65536),
139			FFT_NAME(ff_cos_131072),
140			};
141
142			#endif /* FFT_FIXED_32 */
143
144			static void fft_permute_c(FFTContext s, FFTComplex z);
145			static void fft_calc_c(FFTContext s, FFTComplex z);
146
147		5282064	static int split_radix_permutation(int i, int n, int inverse)
148			{
149			int m;
150	✓✓	5282064	if(n <= 2) return i&1;
151		4505064	m = n >> 1;
152	✓✓	4505064	if(!(i&m)) return split_radix_permutation(i, m, inverse)*2;
153		2252532	m >>= 1;
154	✓✓	2252532	if(inverse == !(i&m)) return split_radix_permutation(i, m, inverse)*4 + 1;
155		1126266	else return split_radix_permutation(i, m, inverse)*4 - 1;
156			}
157
158		10551	av_cold void ff_init_ff_cos_tabs(int index)
159			{
160			#if (!CONFIG_HARDCODED_TABLES) && (!FFT_FIXED_32)
161		10551	ff_thread_once(&cos_tabs_init_once[index].control, cos_tabs_init_once[index].func);
162			#endif
163		10551	}
164
165			static const int avx_tab[] = {
166			0, 4, 1, 5, 8, 12, 9, 13, 2, 6, 3, 7, 10, 14, 11, 15
167			};
168
169		14530	static int is_second_half_of_fft32(int i, int n)
170			{
171	✓✓	14530	if (n <= 32)
172		4090	return i >= 16;
173	✓✓	10440	else if (i < n/2)
174		5220	return is_second_half_of_fft32(i, n/2);
175	✓✓	5220	else if (i < 3*n/4)
176		2610	return is_second_half_of_fft32(i - n/2, n/4);
177			else
178		2610	return is_second_half_of_fft32(i - 3*n/4, n/4);
179			}
180
181		303	static av_cold void fft_perm_avx(FFTContext *s)
182			{
183			int i;
184		303	int n = 1 << s->nbits;
185
186	✓✓	4393	for (i = 0; i < n; i += 16) {
187			int k;
188	✓✓	4090	if (is_second_half_of_fft32(i, n)) {
189	✓✓	23783	for (k = 0; k < 16; k++)
190		44768	s->revtab[-split_radix_permutation(i + k, n, s->inverse) & (n - 1)] =
191		22384	i + avx_tab[k];
192
193			} else {
194	✓✓	45747	for (k = 0; k < 16; k++) {
195		43056	int j = i + k;
196		43056	j = (j & ~7) \| ((j >> 1) & 3) \| ((j << 2) & 4);
197		43056	s->revtab[-split_radix_permutation(i + k, n, s->inverse) & (n - 1)] = j;
198			}
199			}
200			}
201		303	}
202
203		3602	av_cold int ff_fft_init(FFTContext *s, int nbits, int inverse)
204			{
205			int i, j, n;
206
207		3602	s->revtab = NULL;
208		3602	s->revtab32 = NULL;
209
210	✓✗✗✓	3602	if (nbits < 2 \|\| nbits > 17)
211			goto fail;
212		3602	s->nbits = nbits;
213		3602	n = 1 << nbits;
214
215	✓✗	3602	if (nbits <= 16) {
216		3602	s->revtab = av_malloc(n * sizeof(uint16_t));
217	✗✓	3602	if (!s->revtab)
218			goto fail;
219			} else {
220			s->revtab32 = av_malloc(n * sizeof(uint32_t));
221			if (!s->revtab32)
222			goto fail;
223			}
224		3602	s->tmp_buf = av_malloc(n * sizeof(FFTComplex));
225	✗✓	3602	if (!s->tmp_buf)
226			goto fail;
227		3602	s->inverse = inverse;
228		3602	s->fft_permutation = FF_FFT_PERM_DEFAULT;
229
230		3602	s->fft_permute = fft_permute_c;
231		3602	s->fft_calc = fft_calc_c;
232			#if CONFIG_MDCT
233		3602	s->imdct_calc = ff_imdct_calc_c;
234		3602	s->imdct_half = ff_imdct_half_c;
235		3602	s->mdct_calc = ff_mdct_calc_c;
236			#endif
237
238			#if FFT_FIXED_32
239			{
240			static AVOnce control = AV_ONCE_INIT;
241		184	ff_thread_once(&control, fft_lut_init);
242			}
243			#else /* FFT_FIXED_32 */
244			#if FFT_FLOAT
245			if (ARCH_AARCH64) ff_fft_init_aarch64(s);
246			if (ARCH_ARM) ff_fft_init_arm(s);
247			if (ARCH_PPC) ff_fft_init_ppc(s);
248		3378	if (ARCH_X86) ff_fft_init_x86(s);
249		3378	if (CONFIG_MDCT) s->mdct_calcw = s->mdct_calc;
250			if (HAVE_MIPSFPU) ff_fft_init_mips(s);
251			#else
252		40	if (CONFIG_MDCT) s->mdct_calcw = ff_mdct_calcw_c;
253			if (ARCH_ARM) ff_fft_fixed_init_arm(s);
254			#endif
255	✓✓	13790	for(j=4; j<=nbits; j++) {
256		10372	ff_init_ff_cos_tabs(j);
257			}
258			#endif /* FFT_FIXED_32 */
259
260
261	✓✓	3602	if (s->fft_permutation == FF_FFT_PERM_AVX) {
262		303	fft_perm_avx(s);
263			} else {
264			#define PROCESS_FFT_PERM_SWAP_LSBS(num) do {\
265			for(i = 0; i < n; i++) {\
266			int k;\
267			j = i;\
268			j = (j & ~3) \| ((j >> 1) & 1) \| ((j << 1) & 2);\
269			k = -split_radix_permutation(i, n, s->inverse) & (n - 1);\
270			s->revtab##num[k] = j;\
271			} \
272			} while(0);
273
274			#define PROCESS_FFT_PERM_DEFAULT(num) do {\
275			for(i = 0; i < n; i++) {\
276			int k;\
277			j = i;\
278			k = -split_radix_permutation(i, n, s->inverse) & (n - 1);\
279			s->revtab##num[k] = j;\
280			} \
281			} while(0);
282
283			#define SPLIT_RADIX_PERMUTATION(num) do { \
284			if (s->fft_permutation == FF_FFT_PERM_SWAP_LSBS) {\
285			PROCESS_FFT_PERM_SWAP_LSBS(num) \
286			} else {\
287			PROCESS_FFT_PERM_DEFAULT(num) \
288			}\
289			} while(0);
290
291	✓✗	3299	if (s->revtab)
292	✓✓✓✓ ✓✓	714859	SPLIT_RADIX_PERMUTATION()
293	✗✓	3299	if (s->revtab32)
294			SPLIT_RADIX_PERMUTATION(32)
295
296			#undef PROCESS_FFT_PERM_DEFAULT
297			#undef PROCESS_FFT_PERM_SWAP_LSBS
298			#undef SPLIT_RADIX_PERMUTATION
299			}
300
301		3602	return 0;
302			fail:
303			av_freep(&s->revtab);
304			av_freep(&s->revtab32);
305			av_freep(&s->tmp_buf);
306			return -1;
307			}
308
309		48047	static void fft_permute_c(FFTContext s, FFTComplex z)
310			{
311			int j, np;
312		48047	const uint16_t *revtab = s->revtab;
313		48047	const uint32_t *revtab32 = s->revtab32;
314		48047	np = 1 << s->nbits;
315			/* TODO: handle split-radix permute in a more optimal way, probably in-place */
316	✓✗	48047	if (revtab) {
317	✓✓	12913999	for(j=0;j<np;j++) s->tmp_buf[revtab[j]] = z[j];
318			} else
319			for(j=0;j<np;j++) s->tmp_buf[revtab32[j]] = z[j];
320
321		48047	memcpy(z, s->tmp_buf, np * sizeof(FFTComplex));
322		48047	}
323
324		3696	av_cold void ff_fft_end(FFTContext *s)
325			{
326		3696	av_freep(&s->revtab);
327		3696	av_freep(&s->revtab32);
328		3696	av_freep(&s->tmp_buf);
329		3696	}
330
331			#if FFT_FIXED_32
332
333		767510	static void fft_calc_c(FFTContext s, FFTComplex z) {
334
335			int nbits, i, n, num_transforms, offset, step;
336			int n4, n2, n34;
337			unsigned tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8;
338			FFTComplex *tmpz;
339		767510	const int fft_size = (1 << s->nbits);
340			int64_t accu;
341
342		767510	num_transforms = (0x2aab >> (16 - s->nbits)) \| 1;
343
344	✓✓	6540670	for (n=0; n<num_transforms; n++){
345		5773160	offset = ff_fft_offsets_lut[n] << 2;
346		5773160	tmpz = z + offset;
347
348		5773160	tmp1 = tmpz[0].re + (unsigned)tmpz[1].re;
349		5773160	tmp5 = tmpz[2].re + (unsigned)tmpz[3].re;
350		5773160	tmp2 = tmpz[0].im + (unsigned)tmpz[1].im;
351		5773160	tmp6 = tmpz[2].im + (unsigned)tmpz[3].im;
352		5773160	tmp3 = tmpz[0].re - (unsigned)tmpz[1].re;
353		5773160	tmp8 = tmpz[2].im - (unsigned)tmpz[3].im;
354		5773160	tmp4 = tmpz[0].im - (unsigned)tmpz[1].im;
355		5773160	tmp7 = tmpz[2].re - (unsigned)tmpz[3].re;
356
357		5773160	tmpz[0].re = tmp1 + tmp5;
358		5773160	tmpz[2].re = tmp1 - tmp5;
359		5773160	tmpz[0].im = tmp2 + tmp6;
360		5773160	tmpz[2].im = tmp2 - tmp6;
361		5773160	tmpz[1].re = tmp3 + tmp8;
362		5773160	tmpz[3].re = tmp3 - tmp8;
363		5773160	tmpz[1].im = tmp4 - tmp7;
364		5773160	tmpz[3].im = tmp4 + tmp7;
365			}
366
367	✓✓	767510	if (fft_size < 8)
368		1	return;
369
370		767509	num_transforms = (num_transforms >> 1) \| 1;
371
372	✓✓	4020228	for (n=0; n<num_transforms; n++){
373		3252719	offset = ff_fft_offsets_lut[n] << 3;
374		3252719	tmpz = z + offset;
375
376		3252719	tmp1 = tmpz[4].re + (unsigned)tmpz[5].re;
377		3252719	tmp3 = tmpz[6].re + (unsigned)tmpz[7].re;
378		3252719	tmp2 = tmpz[4].im + (unsigned)tmpz[5].im;
379		3252719	tmp4 = tmpz[6].im + (unsigned)tmpz[7].im;
380		3252719	tmp5 = tmp1 + tmp3;
381		3252719	tmp7 = tmp1 - tmp3;
382		3252719	tmp6 = tmp2 + tmp4;
383		3252719	tmp8 = tmp2 - tmp4;
384
385		3252719	tmp1 = tmpz[4].re - (unsigned)tmpz[5].re;
386		3252719	tmp2 = tmpz[4].im - (unsigned)tmpz[5].im;
387		3252719	tmp3 = tmpz[6].re - (unsigned)tmpz[7].re;
388		3252719	tmp4 = tmpz[6].im - (unsigned)tmpz[7].im;
389
390		3252719	tmpz[4].re = tmpz[0].re - tmp5;
391		3252719	tmpz[0].re = tmpz[0].re + tmp5;
392		3252719	tmpz[4].im = tmpz[0].im - tmp6;
393		3252719	tmpz[0].im = tmpz[0].im + tmp6;
394		3252719	tmpz[6].re = tmpz[2].re - tmp8;
395		3252719	tmpz[2].re = tmpz[2].re + tmp8;
396		3252719	tmpz[6].im = tmpz[2].im + tmp7;
397		3252719	tmpz[2].im = tmpz[2].im - tmp7;
398
399		3252719	accu = (int64_t)Q31(M_SQRT1_2)*(int)(tmp1 + tmp2);
400		3252719	tmp5 = (int32_t)((accu + 0x40000000) >> 31);
401		3252719	accu = (int64_t)Q31(M_SQRT1_2)*(int)(tmp3 - tmp4);
402		3252719	tmp7 = (int32_t)((accu + 0x40000000) >> 31);
403		3252719	accu = (int64_t)Q31(M_SQRT1_2)*(int)(tmp2 - tmp1);
404		3252719	tmp6 = (int32_t)((accu + 0x40000000) >> 31);
405		3252719	accu = (int64_t)Q31(M_SQRT1_2)*(int)(tmp3 + tmp4);
406		3252719	tmp8 = (int32_t)((accu + 0x40000000) >> 31);
407		3252719	tmp1 = tmp5 + tmp7;
408		3252719	tmp3 = tmp5 - tmp7;
409		3252719	tmp2 = tmp6 + tmp8;
410		3252719	tmp4 = tmp6 - tmp8;
411
412		3252719	tmpz[5].re = tmpz[1].re - tmp1;
413		3252719	tmpz[1].re = tmpz[1].re + tmp1;
414		3252719	tmpz[5].im = tmpz[1].im - tmp2;
415		3252719	tmpz[1].im = tmpz[1].im + tmp2;
416		3252719	tmpz[7].re = tmpz[3].re - tmp4;
417		3252719	tmpz[3].re = tmpz[3].re + tmp4;
418		3252719	tmpz[7].im = tmpz[3].im + tmp3;
419		3252719	tmpz[3].im = tmpz[3].im - tmp3;
420			}
421
422		767509	step = 1 << ((MAX_LOG2_NFFT-4) - 4);
423		767509	n4 = 4;
424
425	✓✓	2424638	for (nbits=4; nbits<=s->nbits; nbits++){
426		1657129	n2 = 2*n4;
427		1657129	n34 = 3*n4;
428		1657129	num_transforms = (num_transforms >> 1) \| 1;
429
430	✓✓	4159954	for (n=0; n<num_transforms; n++){
431		2502825	const FFTSample *w_re_ptr = ff_w_tab_sr + step;
432		2502825	const FFTSample w_im_ptr = ff_w_tab_sr + MAX_FFT_SIZE/(416) - step;
433		2502825	offset = ff_fft_offsets_lut[n] << nbits;
434		2502825	tmpz = z + offset;
435
436		2502825	tmp5 = tmpz[ n2].re + (unsigned)tmpz[n34].re;
437		2502825	tmp1 = tmpz[ n2].re - (unsigned)tmpz[n34].re;
438		2502825	tmp6 = tmpz[ n2].im + (unsigned)tmpz[n34].im;
439		2502825	tmp2 = tmpz[ n2].im - (unsigned)tmpz[n34].im;
440
441		2502825	tmpz[ n2].re = tmpz[ 0].re - tmp5;
442		2502825	tmpz[ 0].re = tmpz[ 0].re + tmp5;
443		2502825	tmpz[ n2].im = tmpz[ 0].im - tmp6;
444		2502825	tmpz[ 0].im = tmpz[ 0].im + tmp6;
445		2502825	tmpz[n34].re = tmpz[n4].re - tmp2;
446		2502825	tmpz[ n4].re = tmpz[n4].re + tmp2;
447		2502825	tmpz[n34].im = tmpz[n4].im + tmp1;
448		2502825	tmpz[ n4].im = tmpz[n4].im - tmp1;
449
450	✓✓	21313368	for (i=1; i<n4; i++){
451		18810543	FFTSample w_re = w_re_ptr[0];
452		18810543	FFTSample w_im = w_im_ptr[0];
453		18810543	accu = (int64_t)w_re*tmpz[ n2+i].re;
454		18810543	accu += (int64_t)w_im*tmpz[ n2+i].im;
455		18810543	tmp1 = (int32_t)((accu + 0x40000000) >> 31);
456		18810543	accu = (int64_t)w_re*tmpz[ n2+i].im;
457		18810543	accu -= (int64_t)w_im*tmpz[ n2+i].re;
458		18810543	tmp2 = (int32_t)((accu + 0x40000000) >> 31);
459		18810543	accu = (int64_t)w_re*tmpz[n34+i].re;
460		18810543	accu -= (int64_t)w_im*tmpz[n34+i].im;
461		18810543	tmp3 = (int32_t)((accu + 0x40000000) >> 31);
462		18810543	accu = (int64_t)w_re*tmpz[n34+i].im;
463		18810543	accu += (int64_t)w_im*tmpz[n34+i].re;
464		18810543	tmp4 = (int32_t)((accu + 0x40000000) >> 31);
465
466		18810543	tmp5 = tmp1 + tmp3;
467		18810543	tmp1 = tmp1 - tmp3;
468		18810543	tmp6 = tmp2 + tmp4;
469		18810543	tmp2 = tmp2 - tmp4;
470
471		18810543	tmpz[ n2+i].re = tmpz[ i].re - tmp5;
472		18810543	tmpz[ i].re = tmpz[ i].re + tmp5;
473		18810543	tmpz[ n2+i].im = tmpz[ i].im - tmp6;
474		18810543	tmpz[ i].im = tmpz[ i].im + tmp6;
475		18810543	tmpz[n34+i].re = tmpz[n4+i].re - tmp2;
476		18810543	tmpz[ n4+i].re = tmpz[n4+i].re + tmp2;
477		18810543	tmpz[n34+i].im = tmpz[n4+i].im + tmp1;
478		18810543	tmpz[ n4+i].im = tmpz[n4+i].im - tmp1;
479
480		18810543	w_re_ptr += step;
481		18810543	w_im_ptr -= step;
482			}
483			}
484		1657129	step >>= 1;
485		1657129	n4 <<= 1;
486			}
487			}
488
489			#else /* FFT_FIXED_32 */
490
491			#define BUTTERFLIES(a0,a1,a2,a3) {\
492			BF(t3, t5, t5, t1);\
493			BF(a2.re, a0.re, a0.re, t5);\
494			BF(a3.im, a1.im, a1.im, t3);\
495			BF(t4, t6, t2, t6);\
496			BF(a3.re, a1.re, a1.re, t4);\
497			BF(a2.im, a0.im, a0.im, t6);\
498			}
499
500			// force loading all the inputs before storing any.
501			// this is slightly slower for small data, but avoids store->load aliasing
502			// for addresses separated by large powers of 2.
503			#define BUTTERFLIES_BIG(a0,a1,a2,a3) {\
504			FFTSample r0=a0.re, i0=a0.im, r1=a1.re, i1=a1.im;\
505			BF(t3, t5, t5, t1);\
506			BF(a2.re, a0.re, r0, t5);\
507			BF(a3.im, a1.im, i1, t3);\
508			BF(t4, t6, t2, t6);\
509			BF(a3.re, a1.re, r1, t4);\
510			BF(a2.im, a0.im, i0, t6);\
511			}
512
513			#define TRANSFORM(a0,a1,a2,a3,wre,wim) {\
514			CMUL(t1, t2, a2.re, a2.im, wre, -wim);\
515			CMUL(t5, t6, a3.re, a3.im, wre, wim);\
516			BUTTERFLIES(a0,a1,a2,a3)\
517			}
518
519			#define TRANSFORM_ZERO(a0,a1,a2,a3) {\
520			t1 = a2.re;\
521			t2 = a2.im;\
522			t5 = a3.re;\
523			t6 = a3.im;\
524			BUTTERFLIES(a0,a1,a2,a3)\
525			}
526
527			/* z[0...8n-1], w[1...2n-1] */
528			#define PASS(name)\
529			static void name(FFTComplex z, const FFTSample wre, unsigned int n)\
530			{\
531			FFTDouble t1, t2, t3, t4, t5, t6;\
532			int o1 = 2*n;\
533			int o2 = 4*n;\
534			int o3 = 6*n;\
535			const FFTSample *wim = wre+o1;\
536			n--;\
537			\
538			TRANSFORM_ZERO(z[0],z[o1],z[o2],z[o3]);\
539			TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
540			do {\
541			z += 2;\
542			wre += 2;\
543			wim -= 2;\
544			TRANSFORM(z[0],z[o1],z[o2],z[o3],wre[0],wim[0]);\
545			TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
546			} while(--n);\
547			}
548
549	✓✓	25874510	PASS(pass)
550			#if !CONFIG_SMALL
551			#undef BUTTERFLIES
552			#define BUTTERFLIES BUTTERFLIES_BIG
553	✓✓	3066827	PASS(pass_big)
554			#endif
555
556			#define DECL_FFT(n,n2,n4)\
557			static void fft##n(FFTComplex *z)\
558			{\
559			fft##n2(z);\
560			fft##n4(z+n4*2);\
561			fft##n4(z+n4*3);\
562			pass(z,FFT_NAME(ff_cos_##n),n4/2);\
563			}
564
565		19004237	static void fft4(FFTComplex *z)
566			{
567			FFTDouble t1, t2, t3, t4, t5, t6, t7, t8;
568
569		19004237	BF(t3, t1, z[0].re, z[1].re);
570		19004237	BF(t8, t6, z[3].re, z[2].re);
571		19004237	BF(z[2].re, z[0].re, t1, t6);
572		19004237	BF(t4, t2, z[0].im, z[1].im);
573		19004237	BF(t7, t5, z[2].im, z[3].im);
574		19004237	BF(z[3].im, z[1].im, t4, t8);
575		19004237	BF(z[3].re, z[1].re, t3, t7);
576		19004237	BF(z[2].im, z[0].im, t2, t5);
577		19004237	}
578
579		9672333	static void fft8(FFTComplex *z)
580			{
581			FFTDouble t1, t2, t3, t4, t5, t6;
582
583		9672333	fft4(z);
584
585		9672333	BF(t1, z[5].re, z[4].re, -z[5].re);
586		9672333	BF(t2, z[5].im, z[4].im, -z[5].im);
587		9672333	BF(t5, z[7].re, z[6].re, -z[7].re);
588		9672333	BF(t6, z[7].im, z[6].im, -z[7].im);
589
590		9672333	BUTTERFLIES(z[0],z[2],z[4],z[6]);
591		9672333	TRANSFORM(z[1],z[3],z[5],z[7],sqrthalf,sqrthalf);
592		9672333	}
593
594			#if !CONFIG_SMALL
595		4154839	static void fft16(FFTComplex *z)
596			{
597			FFTDouble t1, t2, t3, t4, t5, t6;
598		4154839	FFTSample cos_16_1 = FFT_NAME(ff_cos_16)[1];
599		4154839	FFTSample cos_16_3 = FFT_NAME(ff_cos_16)[3];
600
601		4154839	fft8(z);
602		4154839	fft4(z+8);
603		4154839	fft4(z+12);
604
605		4154839	TRANSFORM_ZERO(z[0],z[4],z[8],z[12]);
606		4154839	TRANSFORM(z[2],z[6],z[10],z[14],sqrthalf,sqrthalf);
607		4154839	TRANSFORM(z[1],z[5],z[9],z[13],cos_16_1,cos_16_3);
608		4154839	TRANSFORM(z[3],z[7],z[11],z[15],cos_16_3,cos_16_1);
609		4154839	}
610			#else
611			DECL_FFT(16,8,4)
612			#endif
613		2604921	DECL_FFT(32,16,8)
614		635379	DECL_FFT(64,32,16)
615		328319	DECL_FFT(128,64,32)
616		114812	DECL_FFT(256,128,64)
617		81399	DECL_FFT(512,256,128)
618			#if !CONFIG_SMALL
619			#define pass pass_big
620			#endif
621		8565	DECL_FFT(1024,512,256)
622		2973	DECL_FFT(2048,1024,512)
623		1010	DECL_FFT(4096,2048,1024)
624		689	DECL_FFT(8192,4096,2048)
625			DECL_FFT(16384,8192,4096)
626			DECL_FFT(32768,16384,8192)
627			DECL_FFT(65536,32768,16384)
628			DECL_FFT(131072,65536,32768)
629
630			static void (* const fft_dispatch[])(FFTComplex*) = {
631			fft4, fft8, fft16, fft32, fft64, fft128, fft256, fft512, fft1024,
632			fft2048, fft4096, fft8192, fft16384, fft32768, fft65536, fft131072
633			};
634
635		3138425	static void fft_calc_c(FFTContext s, FFTComplex z)
636			{
637		3138425	fft_dispatch[s->nbits-2](z);
638		3138425	}
639			#endif /* FFT_FIXED_32 */