mirror of
https://github.com/tildearrow/furnace.git
synced 2024-12-05 18:57:27 +00:00
235 lines
5.4 KiB
C
235 lines
5.4 KiB
C
|
/*
|
||
|
* Copyright (c) 2003, 2007-14 Matteo Frigo
|
||
|
* Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
|
||
|
*
|
||
|
* This program is free software; you can redistribute it and/or modify
|
||
|
* it under the terms of the GNU General Public License as published by
|
||
|
* the Free Software Foundation; either version 2 of the License, or
|
||
|
* (at your option) any later version.
|
||
|
*
|
||
|
* This program 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 General Public License for more details.
|
||
|
*
|
||
|
* You should have received a copy of the GNU General Public License
|
||
|
* along with this program; if not, write to the Free Software
|
||
|
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
|
||
|
*
|
||
|
*/
|
||
|
|
||
|
|
||
|
/* trigonometric functions */
|
||
|
#include "kernel/ifftw.h"
|
||
|
#include <math.h>
|
||
|
|
||
|
#if defined(TRIGREAL_IS_LONG_DOUBLE)
|
||
|
# define COS cosl
|
||
|
# define SIN sinl
|
||
|
# define KTRIG(x) (x##L)
|
||
|
# if defined(HAVE_DECL_SINL) && !HAVE_DECL_SINL
|
||
|
extern long double sinl(long double x);
|
||
|
# endif
|
||
|
# if defined(HAVE_DECL_COSL) && !HAVE_DECL_COSL
|
||
|
extern long double cosl(long double x);
|
||
|
# endif
|
||
|
#elif defined(TRIGREAL_IS_QUAD)
|
||
|
# define COS cosq
|
||
|
# define SIN sinq
|
||
|
# define KTRIG(x) (x##Q)
|
||
|
extern __float128 sinq(__float128 x);
|
||
|
extern __float128 cosq(__float128 x);
|
||
|
#else
|
||
|
# define COS cos
|
||
|
# define SIN sin
|
||
|
# define KTRIG(x) (x)
|
||
|
#endif
|
||
|
|
||
|
static const trigreal K2PI =
|
||
|
KTRIG(6.2831853071795864769252867665590057683943388);
|
||
|
#define by2pi(m, n) ((K2PI * (m)) / (n))
|
||
|
|
||
|
/*
|
||
|
* Improve accuracy by reducing x to range [0..1/8]
|
||
|
* before multiplication by 2 * PI.
|
||
|
*/
|
||
|
|
||
|
static void real_cexp(INT m, INT n, trigreal *out)
|
||
|
{
|
||
|
trigreal theta, c, s, t;
|
||
|
unsigned octant = 0;
|
||
|
INT quarter_n = n;
|
||
|
|
||
|
n += n; n += n;
|
||
|
m += m; m += m;
|
||
|
|
||
|
if (m < 0) m += n;
|
||
|
if (m > n - m) { m = n - m; octant |= 4; }
|
||
|
if (m - quarter_n > 0) { m = m - quarter_n; octant |= 2; }
|
||
|
if (m > quarter_n - m) { m = quarter_n - m; octant |= 1; }
|
||
|
|
||
|
theta = by2pi(m, n);
|
||
|
c = COS(theta); s = SIN(theta);
|
||
|
|
||
|
if (octant & 1) { t = c; c = s; s = t; }
|
||
|
if (octant & 2) { t = c; c = -s; s = t; }
|
||
|
if (octant & 4) { s = -s; }
|
||
|
|
||
|
out[0] = c;
|
||
|
out[1] = s;
|
||
|
}
|
||
|
|
||
|
static INT choose_twshft(INT n)
|
||
|
{
|
||
|
INT log2r = 0;
|
||
|
while (n > 0) {
|
||
|
++log2r;
|
||
|
n /= 4;
|
||
|
}
|
||
|
return log2r;
|
||
|
}
|
||
|
|
||
|
static void cexpl_sqrtn_table(triggen *p, INT m, trigreal *res)
|
||
|
{
|
||
|
m += p->n * (m < 0);
|
||
|
|
||
|
{
|
||
|
INT m0 = m & p->twmsk;
|
||
|
INT m1 = m >> p->twshft;
|
||
|
trigreal wr0 = p->W0[2 * m0];
|
||
|
trigreal wi0 = p->W0[2 * m0 + 1];
|
||
|
trigreal wr1 = p->W1[2 * m1];
|
||
|
trigreal wi1 = p->W1[2 * m1 + 1];
|
||
|
|
||
|
res[0] = wr1 * wr0 - wi1 * wi0;
|
||
|
res[1] = wi1 * wr0 + wr1 * wi0;
|
||
|
}
|
||
|
}
|
||
|
|
||
|
/* multiply (xr, xi) by exp(FFT_SIGN * 2*pi*i*m/n) */
|
||
|
static void rotate_sqrtn_table(triggen *p, INT m, R xr, R xi, R *res)
|
||
|
{
|
||
|
m += p->n * (m < 0);
|
||
|
|
||
|
{
|
||
|
INT m0 = m & p->twmsk;
|
||
|
INT m1 = m >> p->twshft;
|
||
|
trigreal wr0 = p->W0[2 * m0];
|
||
|
trigreal wi0 = p->W0[2 * m0 + 1];
|
||
|
trigreal wr1 = p->W1[2 * m1];
|
||
|
trigreal wi1 = p->W1[2 * m1 + 1];
|
||
|
trigreal wr = wr1 * wr0 - wi1 * wi0;
|
||
|
trigreal wi = wi1 * wr0 + wr1 * wi0;
|
||
|
|
||
|
#if FFT_SIGN == -1
|
||
|
res[0] = xr * wr + xi * wi;
|
||
|
res[1] = xi * wr - xr * wi;
|
||
|
#else
|
||
|
res[0] = xr * wr - xi * wi;
|
||
|
res[1] = xi * wr + xr * wi;
|
||
|
#endif
|
||
|
}
|
||
|
}
|
||
|
|
||
|
static void cexpl_sincos(triggen *p, INT m, trigreal *res)
|
||
|
{
|
||
|
real_cexp(m, p->n, res);
|
||
|
}
|
||
|
|
||
|
static void cexp_zero(triggen *p, INT m, R *res)
|
||
|
{
|
||
|
UNUSED(p); UNUSED(m);
|
||
|
res[0] = 0;
|
||
|
res[1] = 0;
|
||
|
}
|
||
|
|
||
|
static void cexpl_zero(triggen *p, INT m, trigreal *res)
|
||
|
{
|
||
|
UNUSED(p); UNUSED(m);
|
||
|
res[0] = 0;
|
||
|
res[1] = 0;
|
||
|
}
|
||
|
|
||
|
static void cexp_generic(triggen *p, INT m, R *res)
|
||
|
{
|
||
|
trigreal resl[2];
|
||
|
p->cexpl(p, m, resl);
|
||
|
res[0] = (R)resl[0];
|
||
|
res[1] = (R)resl[1];
|
||
|
}
|
||
|
|
||
|
static void rotate_generic(triggen *p, INT m, R xr, R xi, R *res)
|
||
|
{
|
||
|
trigreal w[2];
|
||
|
p->cexpl(p, m, w);
|
||
|
res[0] = xr * w[0] - xi * (FFT_SIGN * w[1]);
|
||
|
res[1] = xi * w[0] + xr * (FFT_SIGN * w[1]);
|
||
|
}
|
||
|
|
||
|
triggen *X(mktriggen)(enum wakefulness wakefulness, INT n)
|
||
|
{
|
||
|
INT i, n0, n1;
|
||
|
triggen *p = (triggen *)MALLOC(sizeof(*p), TWIDDLES);
|
||
|
|
||
|
p->n = n;
|
||
|
p->W0 = p->W1 = 0;
|
||
|
p->cexp = 0;
|
||
|
p->rotate = 0;
|
||
|
|
||
|
switch (wakefulness) {
|
||
|
case SLEEPY:
|
||
|
A(0 /* can't happen */);
|
||
|
break;
|
||
|
|
||
|
case AWAKE_SQRTN_TABLE: {
|
||
|
INT twshft = choose_twshft(n);
|
||
|
|
||
|
p->twshft = twshft;
|
||
|
p->twradix = ((INT)1) << twshft;
|
||
|
p->twmsk = p->twradix - 1;
|
||
|
|
||
|
n0 = p->twradix;
|
||
|
n1 = (n + n0 - 1) / n0;
|
||
|
|
||
|
p->W0 = (trigreal *)MALLOC(n0 * 2 * sizeof(trigreal), TWIDDLES);
|
||
|
p->W1 = (trigreal *)MALLOC(n1 * 2 * sizeof(trigreal), TWIDDLES);
|
||
|
|
||
|
for (i = 0; i < n0; ++i)
|
||
|
real_cexp(i, n, p->W0 + 2 * i);
|
||
|
|
||
|
for (i = 0; i < n1; ++i)
|
||
|
real_cexp(i * p->twradix, n, p->W1 + 2 * i);
|
||
|
|
||
|
p->cexpl = cexpl_sqrtn_table;
|
||
|
p->rotate = rotate_sqrtn_table;
|
||
|
break;
|
||
|
}
|
||
|
|
||
|
case AWAKE_SINCOS:
|
||
|
p->cexpl = cexpl_sincos;
|
||
|
break;
|
||
|
|
||
|
case AWAKE_ZERO:
|
||
|
p->cexp = cexp_zero;
|
||
|
p->cexpl = cexpl_zero;
|
||
|
break;
|
||
|
}
|
||
|
|
||
|
if (!p->cexp) {
|
||
|
if (sizeof(trigreal) == sizeof(R))
|
||
|
p->cexp = (void (*)(triggen *, INT, R *))p->cexpl;
|
||
|
else
|
||
|
p->cexp = cexp_generic;
|
||
|
}
|
||
|
if (!p->rotate)
|
||
|
p->rotate = rotate_generic;
|
||
|
return p;
|
||
|
}
|
||
|
|
||
|
void X(triggen_destroy)(triggen *p)
|
||
|
{
|
||
|
X(ifree0)(p->W0);
|
||
|
X(ifree0)(p->W1);
|
||
|
X(ifree)(p);
|
||
|
}
|