From 14191b3ccd9020928a6e62908534d41565b92ee3 Mon Sep 17 00:00:00 2001
From: Jean-Marc Valin <Jean-Marc.Valin@csiro.au>
Date: Fri, 30 Nov 2007 12:15:49 +1100
Subject: [PATCH] Added pitch analysis. Doesn't crash, but otherwise untested.
---
autogen.sh | 5 -
libcelt/Makefile.am | 4 +-
libcelt/arch.h | 4 +-
libcelt/celt.c | 38 +-
libcelt/celt.h | 4 +-
libcelt/fftwrap.c | 300 +++++++++++
libcelt/fftwrap.h | 58 ++
libcelt/pitch.c | 77 +++
libcelt/pitch.h | 28 +
libcelt/smallft.c | 1260 +++++++++++++++++++++++++++++++++++++++++++
libcelt/smallft.h | 46 ++
libcelt/testcelt.c | 1 +
12 files changed, 1809 insertions(+), 16 deletions(-)
create mode 100644 libcelt/fftwrap.c
create mode 100644 libcelt/fftwrap.h
create mode 100644 libcelt/pitch.c
create mode 100644 libcelt/pitch.h
create mode 100644 libcelt/smallft.c
create mode 100644 libcelt/smallft.h
diff --git a/autogen.sh b/autogen.sh
index 67c73ea48..9abad563e 100755
--- a/autogen.sh
+++ b/autogen.sh
@@ -94,11 +94,6 @@ if test "$DIE" -eq 1; then
exit 1
fi
-if test -z "$*"; then
- echo "I am going to run ./configure with no arguments - if you wish "
- echo "to pass any to it, please specify them on the $0 command line."
-fi
-
echo "Generating configuration files for $package, please wait...."
echo " $ACLOCAL $ACLOCAL_FLAGS"
diff --git a/libcelt/Makefile.am b/libcelt/Makefile.am
index ec58f36c2..d61ed434a 100644
--- a/libcelt/Makefile.am
+++ b/libcelt/Makefile.am
@@ -9,13 +9,13 @@
lib_LTLIBRARIES = libcelt.la
# Sources for compilation in the library
-libcelt_la_SOURCES = celt.c mdct.c
+libcelt_la_SOURCES = celt.c fftwrap.c mdct.c pitch.c smallft.c
#noinst_HEADERS =
libcelt_la_LDFLAGS = -version-info @CELT_LT_CURRENT@:@CELT_LT_REVISION@:@CELT_LT_AGE@
-noinst_HEADERS = arch.h celt.h mdct.h os_support.h
+noinst_HEADERS = arch.h celt.h fftwrap.h mdct.h os_support.h pitch.h smallft.h
noinst_PROGRAMS = testcelt
testcelt_SOURCES = testcelt.c
diff --git a/libcelt/arch.h b/libcelt/arch.h
index ec3c901bc..a398a3963 100644
--- a/libcelt/arch.h
+++ b/libcelt/arch.h
@@ -61,7 +61,7 @@
#else
-#ifndef FLOATING_POINT
+/*#ifndef FLOATING_POINT
#error You now need to define either FIXED_POINT or FLOATING_POINT
#endif
#if defined (ARM4_ASM) || defined(ARM5E_ASM) || defined(BFIN_ASM)
@@ -69,7 +69,7 @@
#endif
#ifdef FIXED_POINT_DEBUG
#error "Don't you think enabling fixed-point is a good thing to do if you want to debug that?"
-#endif
+#endif*/
#endif
diff --git a/libcelt/celt.c b/libcelt/celt.c
index 22c3fd023..75e601e2c 100644
--- a/libcelt/celt.c
+++ b/libcelt/celt.c
@@ -33,9 +33,9 @@
#include "mdct.h"
#include <math.h>
#include "celt.h"
+#include "pitch.h"
-
-#define MAX_PERIOD 2048
+#define MAX_PERIOD 1024
struct CELTState_ {
int frame_size;
@@ -77,13 +77,28 @@ CELTState *celt_encoder_new(int blockSize, int blocksPerFrame)
return st;
}
-void celt_encode(CELTState *st, short *pcm)
+void celt_encoder_destroy(CELTState *st)
+{
+ if (st == NULL)
+ {
+ celt_warning("NULL passed to celt_encoder_destroy");
+ return;
+ }
+ celt_free(st->window);
+ celt_free(st->in_mem);
+ celt_free(st->mdct_overlap);
+ celt_free(st->out_mem);
+ celt_free(st);
+}
+
+int celt_encode(CELTState *st, short *pcm)
{
int i, N, B;
N = st->block_size;
B = st->nb_blocks;
float in[(B+1)*N];
float X[B*N];
+ int pitch_index;
/* FIXME: Add preemphasis */
for (i=0;i<N;i++)
@@ -99,13 +114,19 @@ void celt_encode(CELTState *st, short *pcm)
{
int j;
float x[2*N];
+ float tmp[N];
for (j=0;j<2*N;j++)
x[j] = st->window[j]*in[i*N+j];
- mdct_forward(&st->mdct_lookup, x, X+N*i);
+ mdct_forward(&st->mdct_lookup, x, tmp);
+ /* Interleaving the sub-frames */
+ for (j=0;j<N;j++)
+ X[B*j+i] = tmp[j];
}
+
/* Pitch analysis */
-
+ find_spectral_pitch(in, st->out_mem, MAX_PERIOD, (B+1)*N, &pitch_index, NULL);
+
/* Band normalisation */
/* Pitch prediction */
@@ -120,7 +141,11 @@ void celt_encode(CELTState *st, short *pcm)
{
int j;
float x[2*N];
- mdct_backward(&st->mdct_lookup, X+N*i, x);
+ float tmp[N];
+ /* De-interleaving the sub-frames */
+ for (j=0;j<N;j++)
+ tmp[j] = X[B*j+i];
+ mdct_backward(&st->mdct_lookup, tmp, x);
for (j=0;j<2*N;j++)
x[j] = st->window[j]*x[j];
for (j=0;j<N;j++)
@@ -132,5 +157,6 @@ void celt_encode(CELTState *st, short *pcm)
pcm[i*N+j] = (short)floor(.5+st->out_mem[MAX_PERIOD+(i-B)*N+j]);
}
+ return 0;
}
diff --git a/libcelt/celt.h b/libcelt/celt.h
index 3fe755841..db1b9d4af 100644
--- a/libcelt/celt.h
+++ b/libcelt/celt.h
@@ -39,7 +39,9 @@ typedef struct CELTState_ CELTState;
CELTState *celt_encoder_new(int blockSize, int blocksPerFrame);
-void celt_encode(CELTState *st, short *pcm);
+void celt_encoder_destroy(CELTState *st);
+
+int celt_encode(CELTState *st, short *pcm);
diff --git a/libcelt/fftwrap.c b/libcelt/fftwrap.c
new file mode 100644
index 000000000..31215acc7
--- /dev/null
+++ b/libcelt/fftwrap.c
@@ -0,0 +1,300 @@
+/* Copyright (C) 2005 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
+
+#define USE_SMALLFT
+/*#define USE_KISS_FFT*/
+
+
+#include "arch.h"
+#include "os_support.h"
+
+#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] = 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] = (in[i] + (1<<(shift-1))) >> shift;
+ }
+}
+#endif
+
+#ifdef USE_SMALLFT
+
+#include "smallft.h"
+#include <math.h>
+
+void *spx_fft_init(int size)
+{
+ struct drft_lookup *table;
+ table = celt_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);
+ celt_free(table);
+}
+
+void spx_fft(void *table, float *in, float *out)
+{
+ if (in==out)
+ {
+ int i;
+ celt_warning("FFT should not be done in-place");
+ float scale = 1./((struct drft_lookup *)table)->n;
+ 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)
+ {
+ int i;
+ celt_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_KISS_FFT)
+
+#include "kiss_fftr.h"
+#include "kiss_fft.h"
+
+struct kiss_config {
+ kiss_fftr_cfg forward;
+ kiss_fftr_cfg backward;
+ kiss_fft_cpx *freq_data;
+ int N;
+};
+
+void *spx_fft_init(int size)
+{
+ struct kiss_config *table;
+ table = speex_alloc(sizeof(struct kiss_config));
+ table->freq_data = speex_alloc(sizeof(kiss_fft_cpx)*((size>>1)+1));
+ 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(t->freq_data);
+ speex_free(table);
+}
+
+#ifdef FIXED_POINT
+
+void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
+{
+ int i;
+ int shift;
+ struct kiss_config *t = (struct kiss_config *)table;
+ shift = maximize_range(in, in, 32000, t->N);
+ kiss_fftr(t->forward, in, t->freq_data);
+ out[0] = t->freq_data[0].r;
+ for (i=1;i<t->N>>1;i++)
+ {
+ out[(i<<1)-1] = t->freq_data[i].r;
+ out[(i<<1)] = t->freq_data[i].i;
+ }
+ out[(i<<1)-1] = t->freq_data[i].r;
+ 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_fftr(t->forward, in, t->freq_data);
+ out[0] = scale*t->freq_data[0].r;
+ for (i=1;i<t->N>>1;i++)
+ {
+ out[(i<<1)-1] = scale*t->freq_data[i].r;
+ out[(i<<1)] = scale*t->freq_data[i].i;
+ }
+ out[(i<<1)-1] = scale*t->freq_data[i].r;
+}
+#endif
+
+void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out)
+{
+ int i;
+ struct kiss_config *t = (struct kiss_config *)table;
+ t->freq_data[0].r = in[0];
+ t->freq_data[0].i = 0;
+ for (i=1;i<t->N>>1;i++)
+ {
+ t->freq_data[i].r = in[(i<<1)-1];
+ t->freq_data[i].i = in[(i<<1)];
+ }
+ t->freq_data[i].r = in[(i<<1)-1];
+ t->freq_data[i].i = 0;
+
+ kiss_fftri(t->backward, t->freq_data, out);
+}
+
+
+#else
+
+#error No other FFT implemented
+
+#endif
+
+
+int fixed_point = 1;
+#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
+ spx_word16_t _in[N];
+ spx_word16_t _out[N];
+ 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 (!fixed_point)
+ {
+ struct drft_lookup t;
+ spx_drft_init(&t, ((struct kiss_config *)table)->N);
+ float 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);
+ }
+}
+
+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
+ spx_word16_t _in[N];
+ spx_word16_t _out[N];
+ 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 (!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);
+ }
+}
+
+#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
diff --git a/libcelt/fftwrap.h b/libcelt/fftwrap.h
new file mode 100644
index 000000000..dfaf48944
--- /dev/null
+++ b/libcelt/fftwrap.h
@@ -0,0 +1,58 @@
+/* Copyright (C) 2005 Jean-Marc Valin
+ File: fftwrap.h
+
+ 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.
+
+*/
+
+#ifndef FFTWRAP_H
+#define FFTWRAP_H
+
+#include "arch.h"
+
+/** Compute tables for an FFT */
+void *spx_fft_init(int size);
+
+/** Destroy tables for an FFT */
+void spx_fft_destroy(void *table);
+
+/** Forward (real to half-complex) transform */
+void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out);
+
+/** Backward (half-complex to real) transform */
+void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out);
+
+/** Forward (real to half-complex) transform of float data */
+void spx_fft_float(void *table, float *in, float *out);
+
+/** Backward (half-complex to real) transform of float data */
+void spx_ifft_float(void *table, float *in, float *out);
+
+#endif
diff --git a/libcelt/pitch.c b/libcelt/pitch.c
new file mode 100644
index 000000000..53a6d861c
--- /dev/null
+++ b/libcelt/pitch.c
@@ -0,0 +1,77 @@
+/**
+ @file pitch.c
+ @brief Pitch analysis
+*/
+
+/* Copyright (C) 2005
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ This library 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
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with this library; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#include <stdio.h>
+#include <math.h>
+#include "fftwrap.h"
+
+
+void find_spectral_pitch(float *x, float *y, int lag, int len, int *pitch, float *curve)
+{
+ //FIXME: Yuck!!!
+ static void *fft;
+
+ if (!fft)
+ fft = spx_fft_init(lag);
+
+ float xx[lag];
+ float X[lag];
+ float Y[lag];
+ int i;
+
+ for (i=0;i<lag;i++)
+ xx[i] = 0;
+ for (i=0;i<len;i++)
+ xx[i] = x[i];
+
+ spx_fft(fft, xx, X);
+ spx_fft(fft, y, Y);
+ X[0] = X[0]*Y[0];
+ for (i=1;i<lag/2;i++)
+ {
+ float n = 1.f/(1e1+sqrt((X[2*i-1]*X[2*i-1] + X[2*i ]*X[2*i ])*(Y[2*i-1]*Y[2*i-1] + Y[2*i ]*Y[2*i ])));
+ //n = 1;
+ //n = 1.f/(1+curve[i]);
+
+ float tmp = X[2*i-1];
+ X[2*i-1] = (X[2*i-1]*Y[2*i-1] + X[2*i ]*Y[2*i ])*n;
+ X[2*i ] = (- X[2*i ]*Y[2*i-1] + tmp*Y[2*i ])*n;
+ }
+ X[len-1] = 0;
+ X[0] = X[len-1] = 0;
+ spx_ifft(fft, X, xx);
+
+ float max_corr=-1;
+ //int pitch;
+ *pitch = -1;
+ for (i=0;i<lag-len;i++)
+ {
+ //printf ("%f ", xx[i]);
+ if (xx[i] > max_corr)
+ {
+ *pitch = i;
+ max_corr = xx[i];
+ }
+ }
+ //printf ("%d\n", *pitch);
+}
diff --git a/libcelt/pitch.h b/libcelt/pitch.h
new file mode 100644
index 000000000..9b1bdf90a
--- /dev/null
+++ b/libcelt/pitch.h
@@ -0,0 +1,28 @@
+/**
+ @file pitch.h
+ @brief Pitch analysis
+ */
+
+/* Copyright (C) 2005
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ This library 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
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with this library; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+#ifndef _PITCH_H
+#define _PITCH_H
+
+void find_spectral_pitch(float *x, float *y, int lag, int len, int *pitch, float *curve);
+
+#endif
diff --git a/libcelt/smallft.c b/libcelt/smallft.c
new file mode 100644
index 000000000..35dd1d263
--- /dev/null
+++ b/libcelt/smallft.c
@@ -0,0 +1,1260 @@
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
+ * *
+ * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 *
+ * by the XIPHOPHORUS Company http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: *unnormalized* fft transform
+ last mod: $Id: smallft.c,v 1.19 2003/10/08 05:12:37 jm Exp $
+
+ ********************************************************************/
+
+/* FFT implementation from OggSquish, minus cosine transforms,
+ * minus all but radix 2/4 case. In Vorbis we only need this
+ * cut-down version.
+ *
+ * To do more than just power-of-two sized vectors, see the full
+ * version I wrote for NetLib.
+ *
+ * Note that the packing is a little strange; rather than the FFT r/i
+ * packing following R_0, I_n, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1,
+ * it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the
+ * FORTRAN version
+ */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <math.h>
+#include "smallft.h"
+#include "os_support.h"
+
+static void drfti1(int n, float *wa, int *ifac){
+ static int ntryh[4] = { 4,2,3,5 };
+ static float tpi = 6.28318530717958648f;
+ float arg,argh,argld,fi;
+ int ntry=0,i,j=-1;
+ int k1, l1, l2, ib;
+ int ld, ii, ip, is, nq, nr;
+ int ido, ipm, nfm1;
+ int nl=n;
+ int nf=0;
+
+ L101:
+ j++;
+ if (j < 4)
+ ntry=ntryh[j];
+ else
+ ntry+=2;
+
+ L104:
+ nq=nl/ntry;
+ nr=nl-ntry*nq;
+ if (nr!=0) goto L101;
+
+ nf++;
+ ifac[nf+1]=ntry;
+ nl=nq;
+ if(ntry!=2)goto L107;
+ if(nf==1)goto L107;
+
+ for (i=1;i<nf;i++){
+ ib=nf-i+1;
+ ifac[ib+1]=ifac[ib];
+ }
+ ifac[2] = 2;
+
+ L107:
+ if(nl!=1)goto L104;
+ ifac[0]=n;
+ ifac[1]=nf;
+ argh=tpi/n;
+ is=0;
+ nfm1=nf-1;
+ l1=1;
+
+ if(nfm1==0)return;
+
+ for (k1=0;k1<nfm1;k1++){
+ ip=ifac[k1+2];
+ ld=0;
+ l2=l1*ip;
+ ido=n/l2;
+ ipm=ip-1;
+
+ for (j=0;j<ipm;j++){
+ ld+=l1;
+ i=is;
+ argld=(float)ld*argh;
+ fi=0.f;
+ for (ii=2;ii<ido;ii+=2){
+ fi+=1.f;
+ arg=fi*argld;
+ wa[i++]=cos(arg);
+ wa[i++]=sin(arg);
+ }
+ is+=ido;
+ }
+ l1=l2;
+ }
+}
+
+static void fdrffti(int n, float *wsave, int *ifac){
+
+ if (n == 1) return;
+ drfti1(n, wsave+n, ifac);
+}
+
+static void dradf2(int ido,int l1,float *cc,float *ch,float *wa1){
+ int i,k;
+ float ti2,tr2;
+ int t0,t1,t2,t3,t4,t5,t6;
+
+ t1=0;
+ t0=(t2=l1*ido);
+ t3=ido<<1;
+ for(k=0;k<l1;k++){
+ ch[t1<<1]=cc[t1]+cc[t2];
+ ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];
+ t1+=ido;
+ t2+=ido;
+ }
+
+ if(ido<2)return;
+ if(ido==2)goto L105;
+
+ t1=0;
+ t2=t0;
+ for(k=0;k<l1;k++){
+ t3=t2;
+ t4=(t1<<1)+(ido<<1);
+ t5=t1;
+ t6=t1+t1;
+ for(i=2;i<ido;i+=2){
+ t3+=2;
+ t4-=2;
+ t5+=2;
+ t6+=2;
+ tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
+ ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
+ ch[t6]=cc[t5]+ti2;
+ ch[t4]=ti2-cc[t5];
+ ch[t6-1]=cc[t5-1]+tr2;
+ ch[t4-1]=cc[t5-1]-tr2;
+ }
+ t1+=ido;
+ t2+=ido;
+ }
+
+ if(ido%2==1)return;
+
+ L105:
+ t3=(t2=(t1=ido)-1);
+ t2+=t0;
+ for(k=0;k<l1;k++){
+ ch[t1]=-cc[t2];
+ ch[t1-1]=cc[t3];
+ t1+=ido<<1;
+ t2+=ido;
+ t3+=ido;
+ }
+}
+
+static void dradf4(int ido,int l1,float *cc,float *ch,float *wa1,
+ float *wa2,float *wa3){
+ static float hsqt2 = .70710678118654752f;
+ int i,k,t0,t1,t2,t3,t4,t5,t6;
+ float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
+ t0=l1*ido;
+
+ t1=t0;
+ t4=t1<<1;
+ t2=t1+(t1<<1);
+ t3=0;
+
+ for(k=0;k<l1;k++){
+ tr1=cc[t1]+cc[t2];
+ tr2=cc[t3]+cc[t4];
+
+ ch[t5=t3<<2]=tr1+tr2;
+ ch[(ido<<2)+t5-1]=tr2-tr1;
+ ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];
+ ch[t5]=cc[t2]-cc[t1];
+
+ t1+=ido;
+ t2+=ido;
+ t3+=ido;
+ t4+=ido;
+ }
+
+ if(ido<2)return;
+ if(ido==2)goto L105;
+
+
+ t1=0;
+ for(k=0;k<l1;k++){
+ t2=t1;
+ t4=t1<<2;
+ t5=(t6=ido<<1)+t4;
+ for(i=2;i<ido;i+=2){
+ t3=(t2+=2);
+ t4+=2;
+ t5-=2;
+
+ t3+=t0;
+ cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
+ ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
+ t3+=t0;
+ cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3];
+ ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1];
+ t3+=t0;
+ cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3];
+ ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1];
+
+ tr1=cr2+cr4;
+ tr4=cr4-cr2;
+ ti1=ci2+ci4;
+ ti4=ci2-ci4;
+
+ ti2=cc[t2]+ci3;
+ ti3=cc[t2]-ci3;
+ tr2=cc[t2-1]+cr3;
+ tr3=cc[t2-1]-cr3;
+
+ ch[t4-1]=tr1+tr2;
+ ch[t4]=ti1+ti2;
+
+ ch[t5-1]=tr3-ti4;
+ ch[t5]=tr4-ti3;
+
+ ch[t4+t6-1]=ti4+tr3;
+ ch[t4+t6]=tr4+ti3;
+
+ ch[t5+t6-1]=tr2-tr1;
+ ch[t5+t6]=ti1-ti2;
+ }
+ t1+=ido;
+ }
+ if(ido&1)return;
+
+ L105:
+
+ t2=(t1=t0+ido-1)+(t0<<1);
+ t3=ido<<2;
+ t4=ido;
+ t5=ido<<1;
+ t6=ido;
+
+ for(k=0;k<l1;k++){
+ ti1=-hsqt2*(cc[t1]+cc[t2]);
+ tr1=hsqt2*(cc[t1]-cc[t2]);
+
+ ch[t4-1]=tr1+cc[t6-1];
+ ch[t4+t5-1]=cc[t6-1]-tr1;
+
+ ch[t4]=ti1-cc[t1+t0];
+ ch[t4+t5]=ti1+cc[t1+t0];
+
+ t1+=ido;
+ t2+=ido;
+ t4+=t3;
+ t6+=ido;
+ }
+}
+
+static void dradfg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
+ float *c2,float *ch,float *ch2,float *wa){
+
+ static float tpi=6.283185307179586f;
+ int idij,ipph,i,j,k,l,ic,ik,is;
+ int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
+ float dc2,ai1,ai2,ar1,ar2,ds2;
+ int nbd;
+ float dcp,arg,dsp,ar1h,ar2h;
+ int idp2,ipp2;
+
+ arg=tpi/(float)ip;
+ dcp=cos(arg);
+ dsp=sin(arg);
+ ipph=(ip+1)>>1;
+ ipp2=ip;
+ idp2=ido;
+ nbd=(ido-1)>>1;
+ t0=l1*ido;
+ t10=ip*ido;
+
+ if(ido==1)goto L119;
+ for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];
+
+ t1=0;
+ for(j=1;j<ip;j++){
+ t1+=t0;
+ t2=t1;
+ for(k=0;k<l1;k++){
+ ch[t2]=c1[t2];
+ t2+=ido;
+ }
+ }
+
+ is=-ido;
+ t1=0;
+ if(nbd>l1){
+ for(j=1;j<ip;j++){
+ t1+=t0;
+ is+=ido;
+ t2= -ido+t1;
+ for(k=0;k<l1;k++){
+ idij=is-1;
+ t2+=ido;
+ t3=t2;
+ for(i=2;i<ido;i+=2){
+ idij+=2;
+ t3+=2;
+ ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
+ ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
+ }
+ }
+ }
+ }else{
+
+ for(j=1;j<ip;j++){
+ is+=ido;
+ idij=is-1;
+ t1+=t0;
+ t2=t1;
+ for(i=2;i<ido;i+=2){
+ idij+=2;
+ t2+=2;
+ t3=t2;
+ for(k=0;k<l1;k++){
+ ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
+ ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
+ t3+=ido;
+ }
+ }
+ }
+ }
+
+ t1=0;
+ t2=ipp2*t0;
+ if(nbd<l1){
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(i=2;i<ido;i+=2){
+ t3+=2;
+ t4+=2;
+ t5=t3-ido;
+ t6=t4-ido;
+ for(k=0;k<l1;k++){
+ t5+=ido;
+ t6+=ido;
+ c1[t5-1]=ch[t5-1]+ch[t6-1];
+ c1[t6-1]=ch[t5]-ch[t6];
+ c1[t5]=ch[t5]+ch[t6];
+ c1[t6]=ch[t6-1]-ch[t5-1];
+ }
+ }
+ }
+ }else{
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(k=0;k<l1;k++){
+ t5=t3;
+ t6=t4;
+ for(i=2;i<ido;i+=2){
+ t5+=2;
+ t6+=2;
+ c1[t5-1]=ch[t5-1]+ch[t6-1];
+ c1[t6-1]=ch[t5]-ch[t6];
+ c1[t5]=ch[t5]+ch[t6];
+ c1[t6]=ch[t6-1]-ch[t5-1];
+ }
+ t3+=ido;
+ t4+=ido;
+ }
+ }
+ }
+
+L119:
+ for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
+
+ t1=0;
+ t2=ipp2*idl1;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1-ido;
+ t4=t2-ido;
+ for(k=0;k<l1;k++){
+ t3+=ido;
+ t4+=ido;
+ c1[t3]=ch[t3]+ch[t4];
+ c1[t4]=ch[t4]-ch[t3];
+ }
+ }
+
+ ar1=1.f;
+ ai1=0.f;
+ t1=0;
+ t2=ipp2*idl1;
+ t3=(ip-1)*idl1;
+ for(l=1;l<ipph;l++){
+ t1+=idl1;
+ t2-=idl1;
+ ar1h=dcp*ar1-dsp*ai1;
+ ai1=dcp*ai1+dsp*ar1;
+ ar1=ar1h;
+ t4=t1;
+ t5=t2;
+ t6=t3;
+ t7=idl1;
+
+ for(ik=0;ik<idl1;ik++){
+ ch2[t4++]=c2[ik]+ar1*c2[t7++];
+ ch2[t5++]=ai1*c2[t6++];
+ }
+
+ dc2=ar1;
+ ds2=ai1;
+ ar2=ar1;
+ ai2=ai1;
+
+ t4=idl1;
+ t5=(ipp2-1)*idl1;
+ for(j=2;j<ipph;j++){
+ t4+=idl1;
+ t5-=idl1;
+
+ ar2h=dc2*ar2-ds2*ai2;
+ ai2=dc2*ai2+ds2*ar2;
+ ar2=ar2h;
+
+ t6=t1;
+ t7=t2;
+ t8=t4;
+ t9=t5;
+ for(ik=0;ik<idl1;ik++){
+ ch2[t6++]+=ar2*c2[t8++];
+ ch2[t7++]+=ai2*c2[t9++];
+ }
+ }
+ }
+
+ t1=0;
+ for(j=1;j<ipph;j++){
+ t1+=idl1;
+ t2=t1;
+ for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];
+ }
+
+ if(ido<l1)goto L132;
+
+ t1=0;
+ t2=0;
+ for(k=0;k<l1;k++){
+ t3=t1;
+ t4=t2;
+ for(i=0;i<ido;i++)cc[t4++]=ch[t3++];
+ t1+=ido;
+ t2+=t10;
+ }
+
+ goto L135;
+
+ L132:
+ for(i=0;i<ido;i++){
+ t1=i;
+ t2=i;
+ for(k=0;k<l1;k++){
+ cc[t2]=ch[t1];
+ t1+=ido;
+ t2+=t10;
+ }
+ }
+
+ L135:
+ t1=0;
+ t2=ido<<1;
+ t3=0;
+ t4=ipp2*t0;
+ for(j=1;j<ipph;j++){
+
+ t1+=t2;
+ t3+=t0;
+ t4-=t0;
+
+ t5=t1;
+ t6=t3;
+ t7=t4;
+
+ for(k=0;k<l1;k++){
+ cc[t5-1]=ch[t6];
+ cc[t5]=ch[t7];
+ t5+=t10;
+ t6+=ido;
+ t7+=ido;
+ }
+ }
+
+ if(ido==1)return;
+ if(nbd<l1)goto L141;
+
+ t1=-ido;
+ t3=0;
+ t4=0;
+ t5=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t2;
+ t3+=t2;
+ t4+=t0;
+ t5-=t0;
+ t6=t1;
+ t7=t3;
+ t8=t4;
+ t9=t5;
+ for(k=0;k<l1;k++){
+ for(i=2;i<ido;i+=2){
+ ic=idp2-i;
+ cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];
+ cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];
+ cc[i+t7]=ch[i+t8]+ch[i+t9];
+ cc[ic+t6]=ch[i+t9]-ch[i+t8];
+ }
+ t6+=t10;
+ t7+=t10;
+ t8+=ido;
+ t9+=ido;
+ }
+ }
+ return;
+
+ L141:
+
+ t1=-ido;
+ t3=0;
+ t4=0;
+ t5=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t2;
+ t3+=t2;
+ t4+=t0;
+ t5-=t0;
+ for(i=2;i<ido;i+=2){
+ t6=idp2+t1-i;
+ t7=i+t3;
+ t8=i+t4;
+ t9=i+t5;
+ for(k=0;k<l1;k++){
+ cc[t7-1]=ch[t8-1]+ch[t9-1];
+ cc[t6-1]=ch[t8-1]-ch[t9-1];
+ cc[t7]=ch[t8]+ch[t9];
+ cc[t6]=ch[t9]-ch[t8];
+ t6+=t10;
+ t7+=t10;
+ t8+=ido;
+ t9+=ido;
+ }
+ }
+ }
+}
+
+static void drftf1(int n,float *c,float *ch,float *wa,int *ifac){
+ int i,k1,l1,l2;
+ int na,kh,nf;
+ int ip,iw,ido,idl1,ix2,ix3;
+
+ nf=ifac[1];
+ na=1;
+ l2=n;
+ iw=n;
+
+ for(k1=0;k1<nf;k1++){
+ kh=nf-k1;
+ ip=ifac[kh+1];
+ l1=l2/ip;
+ ido=n/l2;
+ idl1=ido*l1;
+ iw-=(ip-1)*ido;
+ na=1-na;
+
+ if(ip!=4)goto L102;
+
+ ix2=iw+ido;
+ ix3=ix2+ido;
+ if(na!=0)
+ dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
+ else
+ dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
+ goto L110;
+
+ L102:
+ if(ip!=2)goto L104;
+ if(na!=0)goto L103;
+
+ dradf2(ido,l1,c,ch,wa+iw-1);
+ goto L110;
+
+ L103:
+ dradf2(ido,l1,ch,c,wa+iw-1);
+ goto L110;
+
+ L104:
+ if(ido==1)na=1-na;
+ if(na!=0)goto L109;
+
+ dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
+ na=1;
+ goto L110;
+
+ L109:
+ dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
+ na=0;
+
+ L110:
+ l2=l1;
+ }
+
+ if(na==1)return;
+
+ for(i=0;i<n;i++)c[i]=ch[i];
+}
+
+static void dradb2(int ido,int l1,float *cc,float *ch,float *wa1){
+ int i,k,t0,t1,t2,t3,t4,t5,t6;
+ float ti2,tr2;
+
+ t0=l1*ido;
+
+ t1=0;
+ t2=0;
+ t3=(ido<<1)-1;
+ for(k=0;k<l1;k++){
+ ch[t1]=cc[t2]+cc[t3+t2];
+ ch[t1+t0]=cc[t2]-cc[t3+t2];
+ t2=(t1+=ido)<<1;
+ }
+
+ if(ido<2)return;
+ if(ido==2)goto L105;
+
+ t1=0;
+ t2=0;
+ for(k=0;k<l1;k++){
+ t3=t1;
+ t5=(t4=t2)+(ido<<1);
+ t6=t0+t1;
+ for(i=2;i<ido;i+=2){
+ t3+=2;
+ t4+=2;
+ t5-=2;
+ t6+=2;
+ ch[t3-1]=cc[t4-1]+cc[t5-1];
+ tr2=cc[t4-1]-cc[t5-1];
+ ch[t3]=cc[t4]-cc[t5];
+ ti2=cc[t4]+cc[t5];
+ ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2;
+ ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2;
+ }
+ t2=(t1+=ido)<<1;
+ }
+
+ if(ido%2==1)return;
+
+L105:
+ t1=ido-1;
+ t2=ido-1;
+ for(k=0;k<l1;k++){
+ ch[t1]=cc[t2]+cc[t2];
+ ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);
+ t1+=ido;
+ t2+=ido<<1;
+ }
+}
+
+static void dradb3(int ido,int l1,float *cc,float *ch,float *wa1,
+ float *wa2){
+ static float taur = -.5f;
+ static float taui = .8660254037844386f;
+ int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
+ float ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;
+ t0=l1*ido;
+
+ t1=0;
+ t2=t0<<1;
+ t3=ido<<1;
+ t4=ido+(ido<<1);
+ t5=0;
+ for(k=0;k<l1;k++){
+ tr2=cc[t3-1]+cc[t3-1];
+ cr2=cc[t5]+(taur*tr2);
+ ch[t1]=cc[t5]+tr2;
+ ci3=taui*(cc[t3]+cc[t3]);
+ ch[t1+t0]=cr2-ci3;
+ ch[t1+t2]=cr2+ci3;
+ t1+=ido;
+ t3+=t4;
+ t5+=t4;
+ }
+
+ if(ido==1)return;
+
+ t1=0;
+ t3=ido<<1;
+ for(k=0;k<l1;k++){
+ t7=t1+(t1<<1);
+ t6=(t5=t7+t3);
+ t8=t1;
+ t10=(t9=t1+t0)+t0;
+
+ for(i=2;i<ido;i+=2){
+ t5+=2;
+ t6-=2;
+ t7+=2;
+ t8+=2;
+ t9+=2;
+ t10+=2;
+ tr2=cc[t5-1]+cc[t6-1];
+ cr2=cc[t7-1]+(taur*tr2);
+ ch[t8-1]=cc[t7-1]+tr2;
+ ti2=cc[t5]-cc[t6];
+ ci2=cc[t7]+(taur*ti2);
+ ch[t8]=cc[t7]+ti2;
+ cr3=taui*(cc[t5-1]-cc[t6-1]);
+ ci3=taui*(cc[t5]+cc[t6]);
+ dr2=cr2-ci3;
+ dr3=cr2+ci3;
+ di2=ci2+cr3;
+ di3=ci2-cr3;
+ ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2;
+ ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2;
+ ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3;
+ ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3;
+ }
+ t1+=ido;
+ }
+}
+
+static void dradb4(int ido,int l1,float *cc,float *ch,float *wa1,
+ float *wa2,float *wa3){
+ static float sqrt2=1.414213562373095f;
+ int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
+ float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
+ t0=l1*ido;
+
+ t1=0;
+ t2=ido<<2;
+ t3=0;
+ t6=ido<<1;
+ for(k=0;k<l1;k++){
+ t4=t3+t6;
+ t5=t1;
+ tr3=cc[t4-1]+cc[t4-1];
+ tr4=cc[t4]+cc[t4];
+ tr1=cc[t3]-cc[(t4+=t6)-1];
+ tr2=cc[t3]+cc[t4-1];
+ ch[t5]=tr2+tr3;
+ ch[t5+=t0]=tr1-tr4;
+ ch[t5+=t0]=tr2-tr3;
+ ch[t5+=t0]=tr1+tr4;
+ t1+=ido;
+ t3+=t2;
+ }
+
+ if(ido<2)return;
+ if(ido==2)goto L105;
+
+ t1=0;
+ for(k=0;k<l1;k++){
+ t5=(t4=(t3=(t2=t1<<2)+t6))+t6;
+ t7=t1;
+ for(i=2;i<ido;i+=2){
+ t2+=2;
+ t3+=2;
+ t4-=2;
+ t5-=2;
+ t7+=2;
+ ti1=cc[t2]+cc[t5];
+ ti2=cc[t2]-cc[t5];
+ ti3=cc[t3]-cc[t4];
+ tr4=cc[t3]+cc[t4];
+ tr1=cc[t2-1]-cc[t5-1];
+ tr2=cc[t2-1]+cc[t5-1];
+ ti4=cc[t3-1]-cc[t4-1];
+ tr3=cc[t3-1]+cc[t4-1];
+ ch[t7-1]=tr2+tr3;
+ cr3=tr2-tr3;
+ ch[t7]=ti2+ti3;
+ ci3=ti2-ti3;
+ cr2=tr1-tr4;
+ cr4=tr1+tr4;
+ ci2=ti1+ti4;
+ ci4=ti1-ti4;
+
+ ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2;
+ ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2;
+ ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3;
+ ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3;
+ ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4;
+ ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4;
+ }
+ t1+=ido;
+ }
+
+ if(ido%2 == 1)return;
+
+ L105:
+
+ t1=ido;
+ t2=ido<<2;
+ t3=ido-1;
+ t4=ido+(ido<<1);
+ for(k=0;k<l1;k++){
+ t5=t3;
+ ti1=cc[t1]+cc[t4];
+ ti2=cc[t4]-cc[t1];
+ tr1=cc[t1-1]-cc[t4-1];
+ tr2=cc[t1-1]+cc[t4-1];
+ ch[t5]=tr2+tr2;
+ ch[t5+=t0]=sqrt2*(tr1-ti1);
+ ch[t5+=t0]=ti2+ti2;
+ ch[t5+=t0]=-sqrt2*(tr1+ti1);
+
+ t3+=ido;
+ t1+=t2;
+ t4+=t2;
+ }
+}
+
+static void dradbg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
+ float *c2,float *ch,float *ch2,float *wa){
+ static float tpi=6.283185307179586f;
+ int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,
+ t11,t12;
+ float dc2,ai1,ai2,ar1,ar2,ds2;
+ int nbd;
+ float dcp,arg,dsp,ar1h,ar2h;
+ int ipp2;
+
+ t10=ip*ido;
+ t0=l1*ido;
+ arg=tpi/(float)ip;
+ dcp=cos(arg);
+ dsp=sin(arg);
+ nbd=(ido-1)>>1;
+ ipp2=ip;
+ ipph=(ip+1)>>1;
+ if(ido<l1)goto L103;
+
+ t1=0;
+ t2=0;
+ for(k=0;k<l1;k++){
+ t3=t1;
+ t4=t2;
+ for(i=0;i<ido;i++){
+ ch[t3]=cc[t4];
+ t3++;
+ t4++;
+ }
+ t1+=ido;
+ t2+=t10;
+ }
+ goto L106;
+
+ L103:
+ t1=0;
+ for(i=0;i<ido;i++){
+ t2=t1;
+ t3=t1;
+ for(k=0;k<l1;k++){
+ ch[t2]=cc[t3];
+ t2+=ido;
+ t3+=t10;
+ }
+ t1++;
+ }
+
+ L106:
+ t1=0;
+ t2=ipp2*t0;
+ t7=(t5=ido<<1);
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ t6=t5;
+ for(k=0;k<l1;k++){
+ ch[t3]=cc[t6-1]+cc[t6-1];
+ ch[t4]=cc[t6]+cc[t6];
+ t3+=ido;
+ t4+=ido;
+ t6+=t10;
+ }
+ t5+=t7;
+ }
+
+ if (ido == 1)goto L116;
+ if(nbd<l1)goto L112;
+
+ t1=0;
+ t2=ipp2*t0;
+ t7=0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+
+ t7+=(ido<<1);
+ t8=t7;
+ for(k=0;k<l1;k++){
+ t5=t3;
+ t6=t4;
+ t9=t8;
+ t11=t8;
+ for(i=2;i<ido;i+=2){
+ t5+=2;
+ t6+=2;
+ t9+=2;
+ t11-=2;
+ ch[t5-1]=cc[t9-1]+cc[t11-1];
+ ch[t6-1]=cc[t9-1]-cc[t11-1];
+ ch[t5]=cc[t9]-cc[t11];
+ ch[t6]=cc[t9]+cc[t11];
+ }
+ t3+=ido;
+ t4+=ido;
+ t8+=t10;
+ }
+ }
+ goto L116;
+
+ L112:
+ t1=0;
+ t2=ipp2*t0;
+ t7=0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ t7+=(ido<<1);
+ t8=t7;
+ t9=t7;
+ for(i=2;i<ido;i+=2){
+ t3+=2;
+ t4+=2;
+ t8+=2;
+ t9-=2;
+ t5=t3;
+ t6=t4;
+ t11=t8;
+ t12=t9;
+ for(k=0;k<l1;k++){
+ ch[t5-1]=cc[t11-1]+cc[t12-1];
+ ch[t6-1]=cc[t11-1]-cc[t12-1];
+ ch[t5]=cc[t11]-cc[t12];
+ ch[t6]=cc[t11]+cc[t12];
+ t5+=ido;
+ t6+=ido;
+ t11+=t10;
+ t12+=t10;
+ }
+ }
+ }
+
+L116:
+ ar1=1.f;
+ ai1=0.f;
+ t1=0;
+ t9=(t2=ipp2*idl1);
+ t3=(ip-1)*idl1;
+ for(l=1;l<ipph;l++){
+ t1+=idl1;
+ t2-=idl1;
+
+ ar1h=dcp*ar1-dsp*ai1;
+ ai1=dcp*ai1+dsp*ar1;
+ ar1=ar1h;
+ t4=t1;
+ t5=t2;
+ t6=0;
+ t7=idl1;
+ t8=t3;
+ for(ik=0;ik<idl1;ik++){
+ c2[t4++]=ch2[t6++]+ar1*ch2[t7++];
+ c2[t5++]=ai1*ch2[t8++];
+ }
+ dc2=ar1;
+ ds2=ai1;
+ ar2=ar1;
+ ai2=ai1;
+
+ t6=idl1;
+ t7=t9-idl1;
+ for(j=2;j<ipph;j++){
+ t6+=idl1;
+ t7-=idl1;
+ ar2h=dc2*ar2-ds2*ai2;
+ ai2=dc2*ai2+ds2*ar2;
+ ar2=ar2h;
+ t4=t1;
+ t5=t2;
+ t11=t6;
+ t12=t7;
+ for(ik=0;ik<idl1;ik++){
+ c2[t4++]+=ar2*ch2[t11++];
+ c2[t5++]+=ai2*ch2[t12++];
+ }
+ }
+ }
+
+ t1=0;
+ for(j=1;j<ipph;j++){
+ t1+=idl1;
+ t2=t1;
+ for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++];
+ }
+
+ t1=0;
+ t2=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(k=0;k<l1;k++){
+ ch[t3]=c1[t3]-c1[t4];
+ ch[t4]=c1[t3]+c1[t4];
+ t3+=ido;
+ t4+=ido;
+ }
+ }
+
+ if(ido==1)goto L132;
+ if(nbd<l1)goto L128;
+
+ t1=0;
+ t2=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(k=0;k<l1;k++){
+ t5=t3;
+ t6=t4;
+ for(i=2;i<ido;i+=2){
+ t5+=2;
+ t6+=2;
+ ch[t5-1]=c1[t5-1]-c1[t6];
+ ch[t6-1]=c1[t5-1]+c1[t6];
+ ch[t5]=c1[t5]+c1[t6-1];
+ ch[t6]=c1[t5]-c1[t6-1];
+ }
+ t3+=ido;
+ t4+=ido;
+ }
+ }
+ goto L132;
+
+ L128:
+ t1=0;
+ t2=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(i=2;i<ido;i+=2){
+ t3+=2;
+ t4+=2;
+ t5=t3;
+ t6=t4;
+ for(k=0;k<l1;k++){
+ ch[t5-1]=c1[t5-1]-c1[t6];
+ ch[t6-1]=c1[t5-1]+c1[t6];
+ ch[t5]=c1[t5]+c1[t6-1];
+ ch[t6]=c1[t5]-c1[t6-1];
+ t5+=ido;
+ t6+=ido;
+ }
+ }
+ }
+
+L132:
+ if(ido==1)return;
+
+ for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
+
+ t1=0;
+ for(j=1;j<ip;j++){
+ t2=(t1+=t0);
+ for(k=0;k<l1;k++){
+ c1[t2]=ch[t2];
+ t2+=ido;
+ }
+ }
+
+ if(nbd>l1)goto L139;
+
+ is= -ido-1;
+ t1=0;
+ for(j=1;j<ip;j++){
+ is+=ido;
+ t1+=t0;
+ idij=is;
+ t2=t1;
+ for(i=2;i<ido;i+=2){
+ t2+=2;
+ idij+=2;
+ t3=t2;
+ for(k=0;k<l1;k++){
+ c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
+ c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
+ t3+=ido;
+ }
+ }
+ }
+ return;
+
+ L139:
+ is= -ido-1;
+ t1=0;
+ for(j=1;j<ip;j++){
+ is+=ido;
+ t1+=t0;
+ t2=t1;
+ for(k=0;k<l1;k++){
+ idij=is;
+ t3=t2;
+ for(i=2;i<ido;i+=2){
+ idij+=2;
+ t3+=2;
+ c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
+ c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
+ }
+ t2+=ido;
+ }
+ }
+}
+
+static void drftb1(int n, float *c, float *ch, float *wa, int *ifac){
+ int i,k1,l1,l2;
+ int na;
+ int nf,ip,iw,ix2,ix3,ido,idl1;
+
+ nf=ifac[1];
+ na=0;
+ l1=1;
+ iw=1;
+
+ for(k1=0;k1<nf;k1++){
+ ip=ifac[k1 + 2];
+ l2=ip*l1;
+ ido=n/l2;
+ idl1=ido*l1;
+ if(ip!=4)goto L103;
+ ix2=iw+ido;
+ ix3=ix2+ido;
+
+ if(na!=0)
+ dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
+ else
+ dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
+ na=1-na;
+ goto L115;
+
+ L103:
+ if(ip!=2)goto L106;
+
+ if(na!=0)
+ dradb2(ido,l1,ch,c,wa+iw-1);
+ else
+ dradb2(ido,l1,c,ch,wa+iw-1);
+ na=1-na;
+ goto L115;
+
+ L106:
+ if(ip!=3)goto L109;
+
+ ix2=iw+ido;
+ if(na!=0)
+ dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1);
+ else
+ dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1);
+ na=1-na;
+ goto L115;
+
+ L109:
+/* The radix five case can be translated later..... */
+/* if(ip!=5)goto L112;
+
+ ix2=iw+ido;
+ ix3=ix2+ido;
+ ix4=ix3+ido;
+ if(na!=0)
+ dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
+ else
+ dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
+ na=1-na;
+ goto L115;
+
+ L112:*/
+ if(na!=0)
+ dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
+ else
+ dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
+ if(ido==1)na=1-na;
+
+ L115:
+ l1=l2;
+ iw+=(ip-1)*ido;
+ }
+
+ if(na==0)return;
+
+ for(i=0;i<n;i++)c[i]=ch[i];
+}
+
+void spx_drft_forward(struct drft_lookup *l,float *data){
+ if(l->n==1)return;
+ drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
+}
+
+void spx_drft_backward(struct drft_lookup *l,float *data){
+ if (l->n==1)return;
+ drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
+}
+
+void spx_drft_init(struct drft_lookup *l,int n)
+{
+ l->n=n;
+ l->trigcache=(float*)celt_alloc(3*n*sizeof(*l->trigcache));
+ l->splitcache=(int*)celt_alloc(32*sizeof(*l->splitcache));
+ fdrffti(n, l->trigcache, l->splitcache);
+}
+
+void spx_drft_clear(struct drft_lookup *l)
+{
+ if(l)
+ {
+ if(l->trigcache)
+ celt_free(l->trigcache);
+ if(l->splitcache)
+ celt_free(l->splitcache);
+ }
+}
diff --git a/libcelt/smallft.h b/libcelt/smallft.h
new file mode 100644
index 000000000..446e2f65b
--- /dev/null
+++ b/libcelt/smallft.h
@@ -0,0 +1,46 @@
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
+ * *
+ * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 *
+ * by the XIPHOPHORUS Company http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: fft transform
+ last mod: $Id: smallft.h,v 1.3 2003/09/16 18:35:45 jm Exp $
+
+ ********************************************************************/
+/**
+ @file smallft.h
+ @brief Discrete Rotational Fourier Transform (DRFT)
+*/
+
+#ifndef _V_SMFT_H_
+#define _V_SMFT_H_
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/** Discrete Rotational Fourier Transform lookup */
+struct drft_lookup{
+ int n;
+ float *trigcache;
+ int *splitcache;
+};
+
+extern void spx_drft_forward(struct drft_lookup *l,float *data);
+extern void spx_drft_backward(struct drft_lookup *l,float *data);
+extern void spx_drft_init(struct drft_lookup *l,int n);
+extern void spx_drft_clear(struct drft_lookup *l);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/libcelt/testcelt.c b/libcelt/testcelt.c
index 5f9812a3b..010db915b 100644
--- a/libcelt/testcelt.c
+++ b/libcelt/testcelt.c
@@ -56,6 +56,7 @@ int main(int argc, char *argv[])
celt_encode(st, in);
fwrite(in, sizeof(short), FRAME_SIZE, fout);
}
+ celt_encoder_destroy(st);
return 0;
}
--
GitLab