ref: 9714f66d5f75ce400fdb7f38db75695695b125fa
dir: /libcelt/c64_fft.c/
/* (c) Copyright 2008/2009 Xiph.Org Foundation */
/*
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
- Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
- Neither the name of the Xiph.org Foundation nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#include "c64_fft.h"
#include "dsp_fft16x16t.h"
#include "dsp_fft32x32s.h"
#include "dsp_ifft32x32.h"
#ifndef PI
# ifdef M_PI
# define PI M_PI
# else
# define PI 3.14159265358979323846
# endif
#endif
/* ======================================================================== */
/* D2S -- Truncate a 'double' to a 'short', with clamping. */
/* ======================================================================== */
static short d2s(double d)
{
if (d >= 32767.0) return 32767;
if (d <= -32768.0) return -32768;
return (short)d;
}
/* ======================================================================== */
/* D2S -- Truncate a 'double' to a 'int', with clamping. */
/* ======================================================================== */
static int d2i(double d)
{
if (d >= 2147483647.0) return (int)0x7FFFFFFF;
if (d <= -2147483648.0) return (int)0x80000000;
return (int)d;
}
/* ======================================================================== */
/* GEN_TWIDDLE -- Generate twiddle factors for TI's custom FFTs. */
/* */
/* USAGE */
/* This routine is called as follows: */
/* */
/* int gen_twiddle(short *w, int n, double scale) */
/* */
/* short *w Pointer to twiddle-factor array */
/* int n Size of FFT */
/* double scale Scale factor to apply to values. */
/* */
/* The routine will generate the twiddle-factors directly into the */
/* array you specify. The array needs to be approximately 2*N */
/* elements long. (The actual size, which is slightly smaller, is */
/* returned by the function.) */
/* ======================================================================== */
int gen_twiddle16(short *w, int n, double scale)
{
int i, j, k;
for (j = 1, k = 0; j < n >> 2; j = j << 2)
{
for (i = 0; i < n >> 2; i += j << 1)
{
w[k + 11] = d2s(scale * cos(6.0 * PI * (i + j) / n));
w[k + 10] = d2s(scale * sin(6.0 * PI * (i + j) / n));
w[k + 9] = d2s(scale * cos(6.0 * PI * (i ) / n));
w[k + 8] = d2s(scale * sin(6.0 * PI * (i ) / n));
w[k + 7] = d2s(scale * cos(4.0 * PI * (i + j) / n));
w[k + 6] = d2s(scale * sin(4.0 * PI * (i + j) / n));
w[k + 5] = d2s(scale * cos(4.0 * PI * (i ) / n));
w[k + 4] = d2s(scale * sin(4.0 * PI * (i ) / n));
w[k + 3] = d2s(scale * cos(2.0 * PI * (i + j) / n));
w[k + 2] = d2s(scale * sin(2.0 * PI * (i + j) / n));
w[k + 1] = d2s(scale * cos(2.0 * PI * (i ) / n));
w[k + 0] = d2s(scale * sin(2.0 * PI * (i ) / n));
k += 12;
}
}
return k;
}
/* ======================================================================== */
/* GEN_TWIDDLE -- Generate twiddle factors for TI's custom FFTs. */
/* */
/* USAGE */
/* This routine is called as follows: */
/* */
/* int gen_twiddle(int *w, int n, double scale) */
/* */
/* int *w Pointer to twiddle-factor array */
/* int n Size of FFT */
/* double scale Scale factor to apply to values. */
/* */
/* The routine will generate the twiddle-factors directly into the */
/* array you specify. The array needs to be approximately 2*N */
/* elements long. (The actual size, which is slightly smaller, is */
/* returned by the function.) */
/* ======================================================================== */
int gen_twiddle32(int *w, int n, double scale)
{
int i, j, k, s=0, t;
for (j = 1, k = 0; j < n >> 2; j = j << 2, s++)
{
for (i = t=0; i < n >> 2; i += j, t++)
{
w[k + 5] = d2i(scale * cos(6.0 * PI * i / n));
w[k + 4] = d2i(scale * sin(6.0 * PI * i / n));
w[k + 3] = d2i(scale * cos(4.0 * PI * i / n));
w[k + 2] = d2i(scale * sin(4.0 * PI * i / n));
w[k + 1] = d2i(scale * cos(2.0 * PI * i / n));
w[k + 0] = d2i(scale * sin(2.0 * PI * i / n));
k += 6;
}
}
return k;
}
#define NBCACHE 3
static c64_fft_t *cache16[NBCACHE] = {NULL,};
static c64_fft_t *cache32[NBCACHE] = {NULL,};
c64_fft_t *c64_fft16_alloc(int length, int x, int y)
{
c64_fft_t *state;
celt_int16_t *w, *iw;
int i, c;
for (c = 0; c < NBCACHE; c++) {
if (cache16[c] == NULL)
break;
if (cache16[c]->nfft == length)
return cache16[c];
}
state = (c64_fft_t *)celt_alloc(sizeof(c64_fft_t));
state->shift = log(length)/log(2) - ceil(log(length)/log(4)-1);
state->nfft = length;
state->twiddle = celt_alloc(length*2*sizeof(celt_int16_t));
state->itwiddle = celt_alloc(length*2*sizeof(celt_int16_t));
gen_twiddle16((celt_int16_t *)state->twiddle, length, 32767.0);
w = (celt_int16_t *)state->twiddle;
iw = (celt_int16_t *)state->itwiddle;
for (i = 0; i < length; i++) {
iw[2*i+0] = w[2*i+0];
iw[2*i+1] = - w[2*i+1];
}
if (c < NBCACHE)
cache16[c++] = state;
if (c < NBCACHE)
cache16[c] = NULL;
return state;
}
c64_fft_t *c64_fft32_alloc(int length, int x, int y)
{
c64_fft_t *state;
int i, c;
for (c = 0; c < NBCACHE; c++) {
if (cache32[c] == NULL)
break;
if (cache32[c]->nfft == length)
return cache32[c];
}
state = (c64_fft_t *)celt_alloc(sizeof(c64_fft_t));
state->shift = log(length)/log(2) - ceil(log(length)/log(4)-1);
state->nfft = length;
state->twiddle = celt_alloc(length*2*sizeof(celt_int32_t));
state->itwiddle = celt_alloc(length*2*sizeof(celt_int32_t));
// Generate the inverse twiddle first because it does not need scaling
gen_twiddle32(state->itwiddle, length, 2147483647.000000000);
for (i = 0; i < length; i++) {
state->twiddle[2*i+0] = state->itwiddle[2*i+0] >> 1;
state->twiddle[2*i+1] = state->itwiddle[2*i+1] >> 1;
}
if (c < NBCACHE)
cache32[c++] = state;
if (c < NBCACHE)
cache32[c] = NULL;
return state;
}
void c64_fft16_free(c64_fft_t *state)
{
c64_fft32_free(state);
}
void c64_fft32_free(c64_fft_t *state)
{
}
void c64_fft16_inplace(c64_fft_t * restrict state, celt_int16_t *X)
{
int i;
VARDECL(celt_int16_t, cin);
VARDECL(celt_int16_t, cout);
SAVE_STACK;
ALLOC(cin, state->nfft*2, celt_int16_t);
ALLOC(cout, state->nfft*2, celt_int16_t);
for (i = 0; i < state->nfft; i++) {
cin[2*i+0] = X[2*i+0];
cin[2*i+1] = X[2*i+1];
}
DSP_fft16x16t((celt_int16_t *)state->twiddle, state->nfft, cin, cout);
for (i = 0; i < state->nfft; i++) {
X[2*i+0] = cout[2*i+0];
X[2*i+1] = cout[2*i+1];
}
RESTORE_STACK;
}
void c64_fft32(c64_fft_t * restrict state, const celt_int32_t *X, celt_int32_t *Y)
{
int i;
VARDECL(celt_int32_t, cin);
SAVE_STACK;
ALLOC(cin, state->nfft*2, celt_int32_t);
for (i = 0; i < state->nfft; i++) {
cin[2*i+0] = X[2*i+0] >> state->shift;
cin[2*i+1] = X[2*i+1] >> state->shift;
}
DSP_fft32x32s(state->twiddle, state->nfft, cin, Y);
RESTORE_STACK;
}
void c64_ifft16(c64_fft_t * restrict state, const celt_int16_t *X, celt_int16_t *Y)
{
int i;
VARDECL(celt_int16_t, cin);
VARDECL(celt_int16_t, cout);
SAVE_STACK;
ALLOC(cin, state->nfft*2, celt_int16_t);
if ((celt_int32_t)Y & 7)
ALLOC(cout, state->nfft*2, celt_int16_t);
else
cout = Y;
for (i = 0; i < state->nfft; i++) {
// No need to scale for this one but still need to save the input
// because the fft is going to change it!
cin[2*i+0] = X[2*i+0];
cin[2*i+1] = X[2*i+1];
}
DSP_fft16x16t((celt_int16_t *)state->itwiddle, state->nfft, cin, cout);
if ((celt_int32_t)Y & 7)
for (i = 0; i < state->nfft; i++) {
Y[2*i+0] = cout[2*i+0];
Y[2*i+1] = cout[2*i+1];
}
RESTORE_STACK;
}
void c64_ifft32(c64_fft_t * restrict state, const celt_int32_t *X, celt_int32_t *Y)
{
int i;
VARDECL(celt_int32_t, cin);
SAVE_STACK;
ALLOC(cin, state->nfft*2, celt_int32_t);
celt_assert(Y & 7 == 0);
for (i = 0; i < state->nfft; i++) {
// No need to scale for this one but still need to save the input
// because the fft is going to change it!
cin[2*i+0] = X[2*i+0];
cin[2*i+1] = X[2*i+1];
}
DSP_ifft32x32(state->itwiddle, state->nfft, cin, Y);
RESTORE_STACK;
}