mirror of
https://gitlab.freedesktop.org/gstreamer/gstreamer.git
synced 2024-07-03 21:25:54 +00:00
443 lines
11 KiB
C
443 lines
11 KiB
C
/*
|
|
* Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
|
|
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
|
|
*
|
|
* SPDX-License-Identifier: BSD-3-Clause
|
|
* See COPYING file for more information.
|
|
*/
|
|
|
|
|
|
#include "_kiss_fft_guts_f64.h"
|
|
/* The guts header contains all the multiplication and addition macros that are defined for
|
|
fixed or floating point complex numbers. It also delares the kf_ internal functions.
|
|
*/
|
|
|
|
static void
|
|
kf_bfly2 (kiss_fft_f64_cpx * Fout,
|
|
const size_t fstride, const kiss_fft_f64_cfg st, int m)
|
|
{
|
|
kiss_fft_f64_cpx *Fout2;
|
|
kiss_fft_f64_cpx *tw1 = st->twiddles;
|
|
kiss_fft_f64_cpx t;
|
|
Fout2 = Fout + m;
|
|
do {
|
|
C_FIXDIV (*Fout, 2);
|
|
C_FIXDIV (*Fout2, 2);
|
|
|
|
C_MUL (t, *Fout2, *tw1);
|
|
tw1 += fstride;
|
|
C_SUB (*Fout2, *Fout, t);
|
|
C_ADDTO (*Fout, t);
|
|
++Fout2;
|
|
++Fout;
|
|
} while (--m);
|
|
}
|
|
|
|
static void
|
|
kf_bfly4 (kiss_fft_f64_cpx * Fout,
|
|
const size_t fstride, const kiss_fft_f64_cfg st, const size_t m)
|
|
{
|
|
kiss_fft_f64_cpx *tw1, *tw2, *tw3;
|
|
kiss_fft_f64_cpx scratch[6];
|
|
size_t k = m;
|
|
const size_t m2 = 2 * m;
|
|
const size_t m3 = 3 * m;
|
|
|
|
|
|
tw3 = tw2 = tw1 = st->twiddles;
|
|
|
|
do {
|
|
C_FIXDIV (*Fout, 4);
|
|
C_FIXDIV (Fout[m], 4);
|
|
C_FIXDIV (Fout[m2], 4);
|
|
C_FIXDIV (Fout[m3], 4);
|
|
|
|
C_MUL (scratch[0], Fout[m], *tw1);
|
|
C_MUL (scratch[1], Fout[m2], *tw2);
|
|
C_MUL (scratch[2], Fout[m3], *tw3);
|
|
|
|
C_SUB (scratch[5], *Fout, scratch[1]);
|
|
C_ADDTO (*Fout, scratch[1]);
|
|
C_ADD (scratch[3], scratch[0], scratch[2]);
|
|
C_SUB (scratch[4], scratch[0], scratch[2]);
|
|
C_SUB (Fout[m2], *Fout, scratch[3]);
|
|
tw1 += fstride;
|
|
tw2 += fstride * 2;
|
|
tw3 += fstride * 3;
|
|
C_ADDTO (*Fout, scratch[3]);
|
|
|
|
if (st->inverse) {
|
|
Fout[m].r = scratch[5].r - scratch[4].i;
|
|
Fout[m].i = scratch[5].i + scratch[4].r;
|
|
Fout[m3].r = scratch[5].r + scratch[4].i;
|
|
Fout[m3].i = scratch[5].i - scratch[4].r;
|
|
} else {
|
|
Fout[m].r = scratch[5].r + scratch[4].i;
|
|
Fout[m].i = scratch[5].i - scratch[4].r;
|
|
Fout[m3].r = scratch[5].r - scratch[4].i;
|
|
Fout[m3].i = scratch[5].i + scratch[4].r;
|
|
}
|
|
++Fout;
|
|
} while (--k);
|
|
}
|
|
|
|
static void
|
|
kf_bfly3 (kiss_fft_f64_cpx * Fout,
|
|
const size_t fstride, const kiss_fft_f64_cfg st, size_t m)
|
|
{
|
|
size_t k = m;
|
|
const size_t m2 = 2 * m;
|
|
kiss_fft_f64_cpx *tw1, *tw2;
|
|
kiss_fft_f64_cpx scratch[5];
|
|
kiss_fft_f64_cpx epi3;
|
|
epi3 = st->twiddles[fstride * m];
|
|
|
|
tw1 = tw2 = st->twiddles;
|
|
|
|
do {
|
|
C_FIXDIV (*Fout, 3);
|
|
C_FIXDIV (Fout[m], 3);
|
|
C_FIXDIV (Fout[m2], 3);
|
|
|
|
C_MUL (scratch[1], Fout[m], *tw1);
|
|
C_MUL (scratch[2], Fout[m2], *tw2);
|
|
|
|
C_ADD (scratch[3], scratch[1], scratch[2]);
|
|
C_SUB (scratch[0], scratch[1], scratch[2]);
|
|
tw1 += fstride;
|
|
tw2 += fstride * 2;
|
|
|
|
Fout[m].r = Fout->r - HALF_OF (scratch[3].r);
|
|
Fout[m].i = Fout->i - HALF_OF (scratch[3].i);
|
|
|
|
C_MULBYSCALAR (scratch[0], epi3.i);
|
|
|
|
C_ADDTO (*Fout, scratch[3]);
|
|
|
|
Fout[m2].r = Fout[m].r + scratch[0].i;
|
|
Fout[m2].i = Fout[m].i - scratch[0].r;
|
|
|
|
Fout[m].r -= scratch[0].i;
|
|
Fout[m].i += scratch[0].r;
|
|
|
|
++Fout;
|
|
} while (--k);
|
|
}
|
|
|
|
static void
|
|
kf_bfly5 (kiss_fft_f64_cpx * Fout,
|
|
const size_t fstride, const kiss_fft_f64_cfg st, int m)
|
|
{
|
|
kiss_fft_f64_cpx *Fout0, *Fout1, *Fout2, *Fout3, *Fout4;
|
|
int u;
|
|
kiss_fft_f64_cpx scratch[13];
|
|
kiss_fft_f64_cpx *twiddles = st->twiddles;
|
|
kiss_fft_f64_cpx *tw;
|
|
kiss_fft_f64_cpx ya, yb;
|
|
ya = twiddles[fstride * m];
|
|
yb = twiddles[fstride * 2 * m];
|
|
|
|
Fout0 = Fout;
|
|
Fout1 = Fout0 + m;
|
|
Fout2 = Fout0 + 2 * m;
|
|
Fout3 = Fout0 + 3 * m;
|
|
Fout4 = Fout0 + 4 * m;
|
|
|
|
tw = st->twiddles;
|
|
for (u = 0; u < m; ++u) {
|
|
C_FIXDIV (*Fout0, 5);
|
|
C_FIXDIV (*Fout1, 5);
|
|
C_FIXDIV (*Fout2, 5);
|
|
C_FIXDIV (*Fout3, 5);
|
|
C_FIXDIV (*Fout4, 5);
|
|
scratch[0] = *Fout0;
|
|
|
|
C_MUL (scratch[1], *Fout1, tw[u * fstride]);
|
|
C_MUL (scratch[2], *Fout2, tw[2 * u * fstride]);
|
|
C_MUL (scratch[3], *Fout3, tw[3 * u * fstride]);
|
|
C_MUL (scratch[4], *Fout4, tw[4 * u * fstride]);
|
|
|
|
C_ADD (scratch[7], scratch[1], scratch[4]);
|
|
C_SUB (scratch[10], scratch[1], scratch[4]);
|
|
C_ADD (scratch[8], scratch[2], scratch[3]);
|
|
C_SUB (scratch[9], scratch[2], scratch[3]);
|
|
|
|
Fout0->r += scratch[7].r + scratch[8].r;
|
|
Fout0->i += scratch[7].i + scratch[8].i;
|
|
|
|
scratch[5].r =
|
|
scratch[0].r + S_MUL (scratch[7].r, ya.r) + S_MUL (scratch[8].r, yb.r);
|
|
scratch[5].i =
|
|
scratch[0].i + S_MUL (scratch[7].i, ya.r) + S_MUL (scratch[8].i, yb.r);
|
|
|
|
scratch[6].r = S_MUL (scratch[10].i, ya.i) + S_MUL (scratch[9].i, yb.i);
|
|
scratch[6].i = -S_MUL (scratch[10].r, ya.i) - S_MUL (scratch[9].r, yb.i);
|
|
|
|
C_SUB (*Fout1, scratch[5], scratch[6]);
|
|
C_ADD (*Fout4, scratch[5], scratch[6]);
|
|
|
|
scratch[11].r =
|
|
scratch[0].r + S_MUL (scratch[7].r, yb.r) + S_MUL (scratch[8].r, ya.r);
|
|
scratch[11].i =
|
|
scratch[0].i + S_MUL (scratch[7].i, yb.r) + S_MUL (scratch[8].i, ya.r);
|
|
scratch[12].r = -S_MUL (scratch[10].i, yb.i) + S_MUL (scratch[9].i, ya.i);
|
|
scratch[12].i = S_MUL (scratch[10].r, yb.i) - S_MUL (scratch[9].r, ya.i);
|
|
|
|
C_ADD (*Fout2, scratch[11], scratch[12]);
|
|
C_SUB (*Fout3, scratch[11], scratch[12]);
|
|
|
|
++Fout0;
|
|
++Fout1;
|
|
++Fout2;
|
|
++Fout3;
|
|
++Fout4;
|
|
}
|
|
}
|
|
|
|
/* perform the butterfly for one stage of a mixed radix FFT */
|
|
static void
|
|
kf_bfly_generic (kiss_fft_f64_cpx * Fout,
|
|
const size_t fstride, const kiss_fft_f64_cfg st, int m, int p)
|
|
{
|
|
int u, k, q1, q;
|
|
kiss_fft_f64_cpx *twiddles = st->twiddles;
|
|
kiss_fft_f64_cpx t;
|
|
int Norig = st->nfft;
|
|
|
|
kiss_fft_f64_cpx *scratch =
|
|
(kiss_fft_f64_cpx *) KISS_FFT_F64_TMP_ALLOC (sizeof (kiss_fft_f64_cpx) *
|
|
p);
|
|
|
|
for (u = 0; u < m; ++u) {
|
|
k = u;
|
|
for (q1 = 0; q1 < p; ++q1) {
|
|
scratch[q1] = Fout[k];
|
|
C_FIXDIV (scratch[q1], p);
|
|
k += m;
|
|
}
|
|
|
|
k = u;
|
|
for (q1 = 0; q1 < p; ++q1) {
|
|
int twidx = 0;
|
|
Fout[k] = scratch[0];
|
|
for (q = 1; q < p; ++q) {
|
|
twidx += fstride * k;
|
|
if (twidx >= Norig)
|
|
twidx -= Norig;
|
|
C_MUL (t, scratch[q], twiddles[twidx]);
|
|
C_ADDTO (Fout[k], t);
|
|
}
|
|
k += m;
|
|
}
|
|
}
|
|
KISS_FFT_F64_TMP_FREE (scratch);
|
|
}
|
|
|
|
static void
|
|
kf_work (kiss_fft_f64_cpx * Fout,
|
|
const kiss_fft_f64_cpx * f,
|
|
const size_t fstride, int in_stride, int *factors,
|
|
const kiss_fft_f64_cfg st)
|
|
{
|
|
kiss_fft_f64_cpx *Fout_beg = Fout;
|
|
const int p = *factors++; /* the radix */
|
|
const int m = *factors++; /* stage's fft length/p */
|
|
const kiss_fft_f64_cpx *Fout_end = Fout + p * m;
|
|
|
|
#ifdef _OPENMP
|
|
// use openmp extensions at the
|
|
// top-level (not recursive)
|
|
if (fstride == 1 && p <= 5 && m != 1) {
|
|
int k;
|
|
|
|
// execute the p different work units in different threads
|
|
# pragma omp parallel for
|
|
for (k = 0; k < p; ++k)
|
|
kf_work (Fout + k * m, f + fstride * in_stride * k, fstride * p,
|
|
in_stride, factors, st);
|
|
// all threads have joined by this point
|
|
|
|
switch (p) {
|
|
case 2:
|
|
kf_bfly2 (Fout, fstride, st, m);
|
|
break;
|
|
case 3:
|
|
kf_bfly3 (Fout, fstride, st, m);
|
|
break;
|
|
case 4:
|
|
kf_bfly4 (Fout, fstride, st, m);
|
|
break;
|
|
case 5:
|
|
kf_bfly5 (Fout, fstride, st, m);
|
|
break;
|
|
default:
|
|
kf_bfly_generic (Fout, fstride, st, m, p);
|
|
break;
|
|
}
|
|
return;
|
|
}
|
|
#endif
|
|
|
|
if (m == 1) {
|
|
do {
|
|
*Fout = *f;
|
|
f += fstride * in_stride;
|
|
} while (++Fout != Fout_end);
|
|
} else {
|
|
do {
|
|
// recursive call:
|
|
// DFT of size m*p performed by doing
|
|
// p instances of smaller DFTs of size m,
|
|
// each one takes a decimated version of the input
|
|
kf_work (Fout, f, fstride * p, in_stride, factors, st);
|
|
f += fstride * in_stride;
|
|
} while ((Fout += m) != Fout_end);
|
|
}
|
|
|
|
Fout = Fout_beg;
|
|
|
|
// recombine the p smaller DFTs
|
|
switch (p) {
|
|
case 2:
|
|
kf_bfly2 (Fout, fstride, st, m);
|
|
break;
|
|
case 3:
|
|
kf_bfly3 (Fout, fstride, st, m);
|
|
break;
|
|
case 4:
|
|
kf_bfly4 (Fout, fstride, st, m);
|
|
break;
|
|
case 5:
|
|
kf_bfly5 (Fout, fstride, st, m);
|
|
break;
|
|
default:
|
|
kf_bfly_generic (Fout, fstride, st, m, p);
|
|
break;
|
|
}
|
|
}
|
|
|
|
/* facbuf is populated by p1,m1,p2,m2, ...
|
|
where
|
|
p[i] * m[i] = m[i-1]
|
|
m0 = n */
|
|
static void
|
|
kf_factor (int n, int *facbuf)
|
|
{
|
|
int p = 4;
|
|
double floor_sqrt;
|
|
floor_sqrt = floor (sqrt ((double) n));
|
|
|
|
/*factor out powers of 4, powers of 2, then any remaining primes */
|
|
do {
|
|
while (n % p) {
|
|
switch (p) {
|
|
case 4:
|
|
p = 2;
|
|
break;
|
|
case 2:
|
|
p = 3;
|
|
break;
|
|
default:
|
|
p += 2;
|
|
break;
|
|
}
|
|
if (p > floor_sqrt)
|
|
p = n; /* no more factors, skip to end */
|
|
}
|
|
n /= p;
|
|
*facbuf++ = p;
|
|
*facbuf++ = n;
|
|
} while (n > 1);
|
|
}
|
|
|
|
/*
|
|
*
|
|
* User-callable function to allocate all necessary storage space for the fft.
|
|
*
|
|
* The return value is a contiguous block of memory, allocated with malloc. As such,
|
|
* It can be freed with free(), rather than a kiss_fft_f64-specific function.
|
|
* */
|
|
kiss_fft_f64_cfg
|
|
kiss_fft_f64_alloc (int nfft, int inverse_fft, void *mem, size_t *lenmem)
|
|
{
|
|
kiss_fft_f64_cfg st = NULL;
|
|
size_t memneeded = sizeof (struct kiss_fft_f64_state)
|
|
+ sizeof (kiss_fft_f64_cpx) * (nfft - 1); /* twiddle factors */
|
|
|
|
if (lenmem == NULL) {
|
|
st = (kiss_fft_f64_cfg) KISS_FFT_F64_MALLOC (memneeded);
|
|
} else {
|
|
if (mem != NULL && *lenmem >= memneeded)
|
|
st = (kiss_fft_f64_cfg) mem;
|
|
*lenmem = memneeded;
|
|
}
|
|
if (st) {
|
|
int i;
|
|
st->nfft = nfft;
|
|
st->inverse = inverse_fft;
|
|
|
|
for (i = 0; i < nfft; ++i) {
|
|
const double pi =
|
|
3.141592653589793238462643383279502884197169399375105820974944;
|
|
double phase = -2 * pi * i / nfft;
|
|
if (st->inverse)
|
|
phase *= -1;
|
|
kf_cexp (st->twiddles + i, phase);
|
|
}
|
|
|
|
kf_factor (nfft, st->factors);
|
|
}
|
|
return st;
|
|
}
|
|
|
|
|
|
void
|
|
kiss_fft_f64_stride (kiss_fft_f64_cfg st, const kiss_fft_f64_cpx * fin,
|
|
kiss_fft_f64_cpx * fout, int in_stride)
|
|
{
|
|
if (fin == fout) {
|
|
//NOTE: this is not really an in-place FFT algorithm.
|
|
//It just performs an out-of-place FFT into a temp buffer
|
|
kiss_fft_f64_cpx *tmpbuf =
|
|
(kiss_fft_f64_cpx *) KISS_FFT_F64_TMP_ALLOC (sizeof (kiss_fft_f64_cpx) *
|
|
st->nfft);
|
|
kf_work (tmpbuf, fin, 1, in_stride, st->factors, st);
|
|
memcpy (fout, tmpbuf, sizeof (kiss_fft_f64_cpx) * st->nfft);
|
|
KISS_FFT_F64_TMP_FREE (tmpbuf);
|
|
} else {
|
|
kf_work (fout, fin, 1, in_stride, st->factors, st);
|
|
}
|
|
}
|
|
|
|
void
|
|
kiss_fft_f64 (kiss_fft_f64_cfg cfg, const kiss_fft_f64_cpx * fin,
|
|
kiss_fft_f64_cpx * fout)
|
|
{
|
|
kiss_fft_f64_stride (cfg, fin, fout, 1);
|
|
}
|
|
|
|
|
|
void
|
|
kiss_fft_f64_cleanup (void)
|
|
{
|
|
// nothing needed any more
|
|
}
|
|
|
|
int
|
|
kiss_fft_f64_next_fast_size (int n)
|
|
{
|
|
while (1) {
|
|
int m = n;
|
|
while ((m % 2) == 0)
|
|
m /= 2;
|
|
while ((m % 3) == 0)
|
|
m /= 3;
|
|
while ((m % 5) == 0)
|
|
m /= 5;
|
|
if (m <= 1)
|
|
break; /* n is completely factorable by twos, threes, and fives */
|
|
n++;
|
|
}
|
|
return n;
|
|
}
|