ref: bb275c41017e17b9b80e8e060d5450097de0289e
parent: 79b361991e1e0d1e464efc0f2a23f179ce392802
author: Jean-Marc Valin <[email protected]>
date: Tue Aug 11 18:12:50 EDT 2009
Making it possible to use the C64x FFT
--- /dev/null
+++ b/libcelt/c64_fft.c
@@ -1,0 +1,344 @@
+/* (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;
+}
+
+
--- /dev/null
+++ b/libcelt/c64_fft.h
@@ -1,0 +1,58 @@
+/* (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.
+*/
+
+#ifndef _dsp_fft_h_
+#define _dsp_fft_h_
+
+#include "config.h"
+
+#include "arch.h"
+#include "os_support.h"
+#include "mathops.h"
+#include "stack_alloc.h"
+
+typedef struct {+ int nfft;
+ int shift;
+ celt_int32_t *twiddle;
+ celt_int32_t *itwiddle;
+} c64_fft_t;
+
+extern c64_fft_t *c64_fft16_alloc(int length, int x, int y);
+extern void c64_fft16_free(c64_fft_t *state);
+extern void c64_fft16_inplace(c64_fft_t *state, celt_int16_t *X);
+extern void c64_ifft16(c64_fft_t *state, const celt_int16_t *X, celt_int16_t *Y);
+
+extern c64_fft_t *c64_fft32_alloc(int length, int x, int y);
+extern void c64_fft32_free(c64_fft_t *state);
+extern void c64_fft32(c64_fft_t *state, const celt_int32_t *X, celt_int32_t *Y);
+extern void c64_ifft32(c64_fft_t *state, const celt_int32_t *X, celt_int32_t *Y);
+
+#endif
--- a/libcelt/kfft_double.h
+++ b/libcelt/kfft_double.h
@@ -32,7 +32,7 @@
#ifndef KFFT_DOUBLE_H
#define KFFT_DOUBLE_H
-#ifdef ENABLE_TI_DSPLIB
+#ifdef ENABLE_TI_DSPLIB55
#include "dsplib.h"
#include "_kiss_fft_guts.h"
@@ -53,7 +53,18 @@
)
-#else /* ENABLE_TI_DSPLIB */
+#elif defined(ENABLE_TI_DSPLIB64)
+
+#include "kiss_fft.h"
+#include "_kiss_fft_guts.h"
+#include "c64_fft.h"
+
+#define cpx32_fft_alloc(length) (kiss_fft_cfg)(c64_fft32_alloc(length, 0, 0))
+#define cpx32_fft_free(state) c64_fft32_free((c64_fft_t *)state)
+#define cpx32_fft(state, X, Y, nx) c64_fft32 ((c64_fft_t *)state, (const celt_int32_t *)(X), (celt_int32_t *)(Y))
+#define cpx32_ifft(state, X, Y, nx) c64_ifft32((c64_fft_t *)state, (const celt_int32_t *)(X), (celt_int32_t *)(Y))
+
+#else /* ENABLE_TI_DSPLIB55/64 */
#include "kiss_fft.h"
#include "_kiss_fft_guts.h"
--- a/libcelt/kfft_single.h
+++ b/libcelt/kfft_single.h
@@ -32,7 +32,7 @@
#ifndef KFFT_SINGLE_H
#define KFFT_SINGLE_H
-#ifdef ENABLE_TI_DSPLIB
+#ifdef ENABLE_TI_DSPLIB55
#include "dsplib.h"