ref: 2e8a3b20d0a8ab6bc9288cbe8ba7ed5cc2cf8ce6
parent: 42667b0a5f4fde51f211e50b5eb78a1d5c98b1d0
author: Jean-Marc Valin <[email protected]>
date: Mon Feb 25 06:49:38 EST 2008
MDCT conversion, part I.
--- a/libcelt/kiss_fftr.c
+++ b/libcelt/kiss_fftr.c
@@ -84,7 +84,6 @@
/* input buffer timedata is stored row-wise */
int k,ncfft;
kiss_fft_cpx f2k,f1k,tdc,tw;
- celt_word32_t f1kr, f1ki, twr, twi;
ncfft = st->substate->nfft;
--- a/libcelt/mdct.c
+++ b/libcelt/mdct.c
@@ -50,11 +50,14 @@
#include "kiss_fft.h"
#include <math.h>
#include "os_support.h"
+#include "_kiss_fft_guts.h"
#ifndef M_PI
#define M_PI 3.14159263
#endif
+#undef S_MUL
+#define S_MUL(a,b) ((a)*(b))
void mdct_init(mdct_lookup *l,int N)
{int i;
@@ -67,7 +70,6 @@
/* We have enough points that sine isn't necessary */
for (i=0;i<N2;i++)
l->trig[i] = cos(2*M_PI*(i+1./8.)/N);
- l->scale = 1./N4;
}
void mdct_clear(mdct_lookup *l)
@@ -76,36 +78,36 @@
celt_free(l->trig);
}
-void mdct_forward(mdct_lookup *l, celt_sig_t *in, celt_sig_t *out)
+void mdct_forward(mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar *out)
{int i;
int N, N2, N4, N8;
- VARDECL(celt_sig_t *f);
+ VARDECL(kiss_fft_scalar *f);
N = l->n;
N2 = N/2;
N4 = N/4;
N8 = N/8;
- ALLOC(f, N2, celt_sig_t);
+ ALLOC(f, N2, kiss_fft_scalar);
/* Consider the input to be compused of four blocks: [a, b, c, d] */
/* Shuffle, fold, pre-rotate (part 1) */
for(i=0;i<N8;i++)
{- float re, im;
+ kiss_fft_scalar re, im;
/* Real part arranged as -d-cR, Imag part arranged as -b+aR*/
re = -.5*(in[N2+N4+2*i] + in[N2+N4-2*i-1]);
im = -.5*(in[N4+2*i] - in[N4-2*i-1]);
- out[2*i] = re*l->trig[i] - im*l->trig[i+N4];
- out[2*i+1] = im*l->trig[i] + re*l->trig[i+N4];
+ out[2*i] = S_MUL(re,l->trig[i]) - S_MUL(im,l->trig[i+N4]);
+ out[2*i+1] = S_MUL(im,l->trig[i]) + S_MUL(re,l->trig[i+N4]);
}
for(;i<N4;i++)
{- float re, im;
+ kiss_fft_scalar re, im;
/* Real part arranged as a-bR, Imag part arranged as -c-dR */
re = .5*(in[2*i-N4] - in[N2+N4-2*i-1]);
im = -.5*(in[N4+2*i] + in[N+N4-2*i-1]);
- out[2*i] = re*l->trig[i] - im*l->trig[i+N4];
- out[2*i+1] = im*l->trig[i] + re*l->trig[i+N4];
+ out[2*i] = S_MUL(re,l->trig[i]) - S_MUL(im,l->trig[i+N4]);
+ out[2*i+1] = S_MUL(im,l->trig[i]) + S_MUL(re,l->trig[i+N4]);
}
/* N/4 complex FFT, which should normally down-scale by 4/N (but doesn't now) */
@@ -114,28 +116,28 @@
/* Post-rotate and apply the scaling if the FFT doesn't to it itself */
for(i=0;i<N4;i++)
{- out[2*i] = -f[2*i+1]*l->trig[i+N4] + f[2*i] *l->trig[i];
- out[N2-1-2*i] = -f[2*i] *l->trig[i+N4] - f[2*i+1]*l->trig[i];
+ out[2*i] = -S_MUL(f[2*i+1],l->trig[i+N4]) + S_MUL(f[2*i] ,l->trig[i]);
+ out[N2-1-2*i] = -S_MUL(f[2*i] ,l->trig[i+N4]) - S_MUL(f[2*i+1],l->trig[i]);
}
}
-void mdct_backward(mdct_lookup *l, celt_sig_t *in, celt_sig_t *out)
+void mdct_backward(mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar *out)
{int i;
int N, N2, N4, N8;
- VARDECL(celt_sig_t *f);
+ VARDECL(kiss_fft_scalar *f);
N = l->n;
N2 = N/2;
N4 = N/4;
N8 = N/8;
- ALLOC(f, N2, celt_sig_t);
+ ALLOC(f, N2, kiss_fft_scalar);
/* Pre-rotate */
for(i=0;i<N4;i++)
{- out[2*i] = -in[N2-2*i-1] * l->trig[i] - in[2*i]*l->trig[i+N4];
- out[2*i+1] = in[N2-2*i-1] * l->trig[i+N4] - in[2*i]*l->trig[i];
+ out[2*i] = -S_MUL(in[N2-2*i-1], l->trig[i]) - S_MUL(in[2*i],l->trig[i+N4]);
+ out[2*i+1] = S_MUL(in[N2-2*i-1], l->trig[i+N4]) - S_MUL(in[2*i],l->trig[i]);
}
/* Inverse N/4 complex FFT. This one should *not* downscale even in fixed-point */
@@ -144,12 +146,12 @@
/* Post-rotate */
for(i=0;i<N4;i++)
{- float re, im;
+ kiss_fft_scalar re, im;
re = f[2*i];
im = f[2*i+1];
/* We'd scale up by 2 here, but instead it's done when mixing the windows */
- f[2*i] = re*l->trig[i] + im*l->trig[i+N4];
- f[2*i+1] = im*l->trig[i] - re*l->trig[i+N4];
+ f[2*i] = S_MUL(re,l->trig[i]) + S_MUL(im,l->trig[i+N4]);
+ f[2*i+1] = S_MUL(im,l->trig[i]) - S_MUL(re,l->trig[i+N4]);
}
/* De-shuffle the components for the middle of the window only */
for(i = 0; i < N4; i++)
--- a/libcelt/mdct.h
+++ b/libcelt/mdct.h
@@ -48,7 +48,6 @@
int n;
kiss_fft_cfg kfft;
float *trig;
- float scale;
} mdct_lookup;
void mdct_init(mdct_lookup *l,int N);
@@ -55,8 +54,8 @@
void mdct_clear(mdct_lookup *l);
/** Compute a forward MDCT and scale by 2/N */
-void mdct_forward(mdct_lookup *l, celt_sig_t *in, celt_sig_t *out);
+void mdct_forward(mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar *out);
/** Compute a backward MDCT (no scaling) */
-void mdct_backward(mdct_lookup *l, celt_sig_t *in, celt_sig_t *out);
+void mdct_backward(mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar *out);
--- a/tests/mdct-test.c
+++ b/tests/mdct-test.c
@@ -6,7 +6,7 @@
#include "mdct.h"
int ret = 0;
-void check(float * in,float * out,int nfft,int isinverse)
+void check(kiss_fft_scalar * in,kiss_fft_scalar * out,int nfft,int isinverse)
{int bin,k;
double errpow=0,sigpow=0;
@@ -36,7 +36,7 @@
}
}
-void check_inv(float * in,float * out,int nfft,int isinverse)
+void check_inv(kiss_fft_scalar * in,kiss_fft_scalar * out,int nfft,int isinverse)
{int bin,k;
double errpow=0,sigpow=0;
@@ -70,10 +70,10 @@
void test1d(int nfft,int isinverse)
{mdct_lookup cfg;
- size_t buflen = sizeof(float)*nfft;
+ size_t buflen = sizeof(kiss_fft_scalar)*nfft;
- float * in = (float*)malloc(buflen);
- float * out= (float*)malloc(buflen);
+ kiss_fft_scalar * in = (kiss_fft_scalar*)malloc(buflen);
+ kiss_fft_scalar * out= (kiss_fft_scalar*)malloc(buflen);
int k;
mdct_init(&cfg, nfft);
@@ -83,8 +83,7 @@
#ifdef DOUBLE_PRECISION
for (k=0;k<nfft;++k) {- in[k].r *= 65536;
- in[k].i *= 65536;
+ in[k] *= 32768;
}
#endif