summaryrefslogtreecommitdiff
path: root/libspeex/fftwrap.c
diff options
context:
space:
mode:
Diffstat (limited to 'libspeex/fftwrap.c')
-rw-r--r--libspeex/fftwrap.c397
1 files changed, 397 insertions, 0 deletions
diff --git a/libspeex/fftwrap.c b/libspeex/fftwrap.c
new file mode 100644
index 0000000..4f37e1b
--- /dev/null
+++ b/libspeex/fftwrap.c
@@ -0,0 +1,397 @@
+/* Copyright (C) 2005-2006 Jean-Marc Valin
+ File: fftwrap.c
+
+ Wrapper for various FFTs
+
+ 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.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "arch.h"
+#include "os_support.h"
+
+#define MAX_FFT_SIZE 2048
+
+#ifdef FIXED_POINT
+static int maximize_range(spx_word16_t *in, spx_word16_t *out, spx_word16_t bound, int len)
+{
+ int i, shift;
+ spx_word16_t max_val = 0;
+ for (i=0;i<len;i++)
+ {
+ if (in[i]>max_val)
+ max_val = in[i];
+ if (-in[i]>max_val)
+ max_val = -in[i];
+ }
+ shift=0;
+ while (max_val <= (bound>>1) && max_val != 0)
+ {
+ max_val <<= 1;
+ shift++;
+ }
+ for (i=0;i<len;i++)
+ {
+ out[i] = SHL16(in[i], shift);
+ }
+ return shift;
+}
+
+static void renorm_range(spx_word16_t *in, spx_word16_t *out, int shift, int len)
+{
+ int i;
+ for (i=0;i<len;i++)
+ {
+ out[i] = PSHR16(in[i], shift);
+ }
+}
+#endif
+
+#ifdef USE_SMALLFT
+
+#include "smallft.h"
+#include <math.h>
+
+void *spx_fft_init(int size)
+{
+ struct drft_lookup *table;
+ table = speex_alloc(sizeof(struct drft_lookup));
+ spx_drft_init((struct drft_lookup *)table, size);
+ return (void*)table;
+}
+
+void spx_fft_destroy(void *table)
+{
+ spx_drft_clear(table);
+ speex_free(table);
+}
+
+void spx_fft(void *table, float *in, float *out)
+{
+ if (in==out)
+ {
+ int i;
+ float scale = 1./((struct drft_lookup *)table)->n;
+ speex_warning("FFT should not be done in-place");
+ for (i=0;i<((struct drft_lookup *)table)->n;i++)
+ out[i] = scale*in[i];
+ } else {
+ int i;
+ float scale = 1./((struct drft_lookup *)table)->n;
+ for (i=0;i<((struct drft_lookup *)table)->n;i++)
+ out[i] = scale*in[i];
+ }
+ spx_drft_forward((struct drft_lookup *)table, out);
+}
+
+void spx_ifft(void *table, float *in, float *out)
+{
+ if (in==out)
+ {
+ speex_warning("FFT should not be done in-place");
+ } else {
+ int i;
+ for (i=0;i<((struct drft_lookup *)table)->n;i++)
+ out[i] = in[i];
+ }
+ spx_drft_backward((struct drft_lookup *)table, out);
+}
+
+#elif defined(USE_INTEL_MKL)
+#include <mkl.h>
+
+struct mkl_config {
+ DFTI_DESCRIPTOR_HANDLE desc;
+ int N;
+};
+
+void *spx_fft_init(int size)
+{
+ struct mkl_config *table = (struct mkl_config *) speex_alloc(sizeof(struct mkl_config));
+ table->N = size;
+ DftiCreateDescriptor(&table->desc, DFTI_SINGLE, DFTI_REAL, 1, size);
+ DftiSetValue(table->desc, DFTI_PACKED_FORMAT, DFTI_PACK_FORMAT);
+ DftiSetValue(table->desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
+ DftiSetValue(table->desc, DFTI_FORWARD_SCALE, 1.0f / size);
+ DftiCommitDescriptor(table->desc);
+ return table;
+}
+
+void spx_fft_destroy(void *table)
+{
+ struct mkl_config *t = (struct mkl_config *) table;
+ DftiFreeDescriptor(t->desc);
+ speex_free(table);
+}
+
+void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
+{
+ struct mkl_config *t = (struct mkl_config *) table;
+ DftiComputeForward(t->desc, in, out);
+}
+
+void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out)
+{
+ struct mkl_config *t = (struct mkl_config *) table;
+ DftiComputeBackward(t->desc, in, out);
+}
+
+#elif defined(USE_GPL_FFTW3)
+
+#include <fftw3.h>
+
+struct fftw_config {
+ float *in;
+ float *out;
+ fftwf_plan fft;
+ fftwf_plan ifft;
+ int N;
+};
+
+void *spx_fft_init(int size)
+{
+ struct fftw_config *table = (struct fftw_config *) speex_alloc(sizeof(struct fftw_config));
+ table->in = fftwf_malloc(sizeof(float) * (size+2));
+ table->out = fftwf_malloc(sizeof(float) * (size+2));
+
+ table->fft = fftwf_plan_dft_r2c_1d(size, table->in, (fftwf_complex *) table->out, FFTW_PATIENT);
+ table->ifft = fftwf_plan_dft_c2r_1d(size, (fftwf_complex *) table->in, table->out, FFTW_PATIENT);
+
+ table->N = size;
+ return table;
+}
+
+void spx_fft_destroy(void *table)
+{
+ struct fftw_config *t = (struct fftw_config *) table;
+ fftwf_destroy_plan(t->fft);
+ fftwf_destroy_plan(t->ifft);
+ fftwf_free(t->in);
+ fftwf_free(t->out);
+ speex_free(table);
+}
+
+
+void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
+{
+ int i;
+ struct fftw_config *t = (struct fftw_config *) table;
+ const int N = t->N;
+ float *iptr = t->in;
+ float *optr = t->out;
+ const float m = 1.0 / N;
+ for(i=0;i<N;++i)
+ iptr[i]=in[i] * m;
+
+ fftwf_execute(t->fft);
+
+ out[0] = optr[0];
+ for(i=1;i<N;++i)
+ out[i] = optr[i+1];
+}
+
+void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out)
+{
+ int i;
+ struct fftw_config *t = (struct fftw_config *) table;
+ const int N = t->N;
+ float *iptr = t->in;
+ float *optr = t->out;
+
+ iptr[0] = in[0];
+ iptr[1] = 0.0f;
+ for(i=1;i<N;++i)
+ iptr[i+1] = in[i];
+ iptr[N+1] = 0.0f;
+
+ fftwf_execute(t->ifft);
+
+ for(i=0;i<N;++i)
+ out[i] = optr[i];
+}
+
+#elif defined(USE_KISS_FFT)
+
+#include "kiss_fftr.h"
+#include "kiss_fft.h"
+
+struct kiss_config {
+ kiss_fftr_cfg forward;
+ kiss_fftr_cfg backward;
+ int N;
+};
+
+void *spx_fft_init(int size)
+{
+ struct kiss_config *table;
+ table = (struct kiss_config*)speex_alloc(sizeof(struct kiss_config));
+ table->forward = kiss_fftr_alloc(size,0,NULL,NULL);
+ table->backward = kiss_fftr_alloc(size,1,NULL,NULL);
+ table->N = size;
+ return table;
+}
+
+void spx_fft_destroy(void *table)
+{
+ struct kiss_config *t = (struct kiss_config *)table;
+ kiss_fftr_free(t->forward);
+ kiss_fftr_free(t->backward);
+ speex_free(table);
+}
+
+#ifdef FIXED_POINT
+
+void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
+{
+ int shift;
+ struct kiss_config *t = (struct kiss_config *)table;
+ shift = maximize_range(in, in, 32000, t->N);
+ kiss_fftr2(t->forward, in, out);
+ renorm_range(in, in, shift, t->N);
+ renorm_range(out, out, shift, t->N);
+}
+
+#else
+
+void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
+{
+ int i;
+ float scale;
+ struct kiss_config *t = (struct kiss_config *)table;
+ scale = 1./t->N;
+ kiss_fftr2(t->forward, in, out);
+ for (i=0;i<t->N;i++)
+ out[i] *= scale;
+}
+#endif
+
+void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out)
+{
+ struct kiss_config *t = (struct kiss_config *)table;
+ kiss_fftri2(t->backward, in, out);
+}
+
+
+#else
+
+#error No other FFT implemented
+
+#endif
+
+
+#ifdef FIXED_POINT
+/*#include "smallft.h"*/
+
+
+void spx_fft_float(void *table, float *in, float *out)
+{
+ int i;
+#ifdef USE_SMALLFT
+ int N = ((struct drft_lookup *)table)->n;
+#elif defined(USE_KISS_FFT)
+ int N = ((struct kiss_config *)table)->N;
+#else
+#endif
+#ifdef VAR_ARRAYS
+ spx_word16_t _in[N];
+ spx_word16_t _out[N];
+#else
+ spx_word16_t _in[MAX_FFT_SIZE];
+ spx_word16_t _out[MAX_FFT_SIZE];
+#endif
+ for (i=0;i<N;i++)
+ _in[i] = (int)floor(.5+in[i]);
+ spx_fft(table, _in, _out);
+ for (i=0;i<N;i++)
+ out[i] = _out[i];
+#if 0
+ if (!fixed_point)
+ {
+ float scale;
+ struct drft_lookup t;
+ spx_drft_init(&t, ((struct kiss_config *)table)->N);
+ scale = 1./((struct kiss_config *)table)->N;
+ for (i=0;i<((struct kiss_config *)table)->N;i++)
+ out[i] = scale*in[i];
+ spx_drft_forward(&t, out);
+ spx_drft_clear(&t);
+ }
+#endif
+}
+
+void spx_ifft_float(void *table, float *in, float *out)
+{
+ int i;
+#ifdef USE_SMALLFT
+ int N = ((struct drft_lookup *)table)->n;
+#elif defined(USE_KISS_FFT)
+ int N = ((struct kiss_config *)table)->N;
+#else
+#endif
+#ifdef VAR_ARRAYS
+ spx_word16_t _in[N];
+ spx_word16_t _out[N];
+#else
+ spx_word16_t _in[MAX_FFT_SIZE];
+ spx_word16_t _out[MAX_FFT_SIZE];
+#endif
+ for (i=0;i<N;i++)
+ _in[i] = (int)floor(.5+in[i]);
+ spx_ifft(table, _in, _out);
+ for (i=0;i<N;i++)
+ out[i] = _out[i];
+#if 0
+ if (!fixed_point)
+ {
+ int i;
+ struct drft_lookup t;
+ spx_drft_init(&t, ((struct kiss_config *)table)->N);
+ for (i=0;i<((struct kiss_config *)table)->N;i++)
+ out[i] = in[i];
+ spx_drft_backward(&t, out);
+ spx_drft_clear(&t);
+ }
+#endif
+}
+
+#else
+
+void spx_fft_float(void *table, float *in, float *out)
+{
+ spx_fft(table, in, out);
+}
+void spx_ifft_float(void *table, float *in, float *out)
+{
+ spx_ifft(table, in, out);
+}
+
+#endif
View Lorry Controller Status