fixed-point: changed find_spectral_pitch() to use single-precision (16-bit) FFT.
authorJean-Marc Valin <Jean-Marc.Valin@csiro.au>
Wed, 5 Mar 2008 06:20:30 +0000 (17:20 +1100)
committerJean-Marc Valin <Jean-Marc.Valin@csiro.au>
Wed, 5 Mar 2008 06:20:30 +0000 (17:20 +1100)
This involved adding kfft_single.[ch] that redefines kiss_fft a second time
with a different prefix. All this is still a bit of a mess now. The mask
had to be converted to 16-bit input, but we're still using floats to apply it.

12 files changed:
libcelt/Makefile.am
libcelt/_kiss_fft_guts.h
libcelt/celt.c
libcelt/kfft_single.c [new file with mode: 0644]
libcelt/kfft_single.h [new file with mode: 0644]
libcelt/kiss_fft.c
libcelt/kiss_fft.h
libcelt/kiss_fftr.c
libcelt/pitch.c
libcelt/pitch.h
libcelt/psy.c
libcelt/psy.h

index 795c8fb..d84b68b 100644 (file)
@@ -14,8 +14,8 @@ noinst_SCRIPTS = match-test.sh
 lib_LTLIBRARIES = libcelt.la
 
 # Sources for compilation in the library
-libcelt_la_SOURCES = bands.c celt.c cwrs.c \
-       ecintrin.h entcode.c entdec.c entenc.c header.c kiss_fft.c \
+libcelt_la_SOURCES = bands.c celt.c cwrs.c ecintrin.h \
+       entcode.c entdec.c entenc.c header.c kfft_single.c kiss_fft.c \
        kiss_fftr.c laplace.c mdct.c modes.c pitch.c psy.c quant_bands.c \
        quant_pitch.c rangedec.c rangeenc.c rate.c vq.c
 
@@ -25,8 +25,8 @@ libcelt_la_LDFLAGS = -version-info @CELT_LT_CURRENT@:@CELT_LT_REVISION@:@CELT_LT
 
 noinst_HEADERS = _kiss_fft_guts.h arch.h bands.h \
        cwrs.h ecintrin.h entcode.h entdec.h entenc.h fixed_generic.h \
-       kiss_fft.h kiss_fftr.h laplace.h mdct.h mfrngcod.h mathops.h modes.h \
-       os_support.h pgain_table.h pitch.h psy.h \
+       kfft_single.h kiss_fft.h kiss_fftr.h laplace.h mdct.h mfrngcod.h \
+       mathops.h modes.h os_support.h pgain_table.h pitch.h psy.h \
        quant_bands.h quant_pitch.h rate.h stack_alloc.h vq.h
 
 noinst_PROGRAMS = testcelt
index ddbe16c..8cc6cba 100644 (file)
@@ -12,6 +12,9 @@ Redistribution and use in source and binary forms, with or without modification,
 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 COPYRIGHT OWNER 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 KISS_FFT_GUTS_H
+#define KISS_FFT_GUTS_H
+
 #define MIN(a,b) ((a)<(b) ? (a):(b))
 #define MAX(a,b) ((a)>(b) ? (a):(b))
 
@@ -229,3 +232,5 @@ struct kiss_fft_state{
 /* a debugging function */
 #define pcpx(c)\
     fprintf(stderr,"%g + %gi\n",(double)((c)->r),(double)((c)->i) )
+
+#endif /* KISS_FFT_GUTS_H */
index e1c015f..7cecb5b 100644 (file)
@@ -103,7 +103,7 @@ CELTEncoder *celt_encoder_create(const CELTMode *mode)
    ec_byte_writeinit(&st->buf);
    ec_enc_init(&st->enc,&st->buf);
 
-   st->fft = kiss_fftr_alloc(MAX_PERIOD, 0, 0);
+   st->fft = pitch_state_alloc(MAX_PERIOD);
    psydecay_init(&st->psy, MAX_PERIOD/2, st->mode->Fs);
    
    st->in_mem = celt_alloc(N*C*sizeof(celt_sig_t));
@@ -130,7 +130,7 @@ void celt_encoder_destroy(CELTEncoder *st)
 
    ec_byte_writeclear(&st->buf);
 
-   kiss_fft_free(st->fft);
+   pitch_state_free(st->fft);
    psydecay_clear(&st->psy);
 
    celt_free(st->in_mem);
diff --git a/libcelt/kfft_single.c b/libcelt/kfft_single.c
new file mode 100644 (file)
index 0000000..9416d45
--- /dev/null
@@ -0,0 +1,13 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef FIXED_POINT
+
+#include "kfft_single.h"
+
+#define SKIP_CONFIG_H
+#include "kiss_fft.c"
+#include "kiss_fftr.c"
+
+#endif
diff --git a/libcelt/kfft_single.h b/libcelt/kfft_single.h
new file mode 100644 (file)
index 0000000..7bc081b
--- /dev/null
@@ -0,0 +1,51 @@
+/* (C) 2008 Jean-Marc Valin, CSIRO
+*/
+/*
+   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 KFFT_SINGLE_H
+#define KFFT_SINGLE_H
+
+#ifdef FIXED_POINT
+
+#ifdef DOUBLE_PRECISION
+#undef DOUBLE_PRECISION
+#endif
+
+#ifdef MIXED_PRECISION
+#undef MIXED_PRECISION
+#endif
+
+#endif /* FIXED_POINT */
+
+#include "kiss_fft.h"
+#include "kiss_fftr.h"
+#include "_kiss_fft_guts.h"
+
+#endif /* KFFT_SINGLE_H */
index d8b7697..b072044 100644 (file)
@@ -15,9 +15,10 @@ Redistribution and use in source and binary forms, with or without modification,
 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 COPYRIGHT OWNER 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"
+#ifndef SKIP_CONFIG_H
+#  ifdef HAVE_CONFIG_H
+#    include "config.h"
+#  endif
 #endif
 
 #include "_kiss_fft_guts.h"
index 568ec6f..f71a839 100644 (file)
@@ -36,22 +36,24 @@ extern "C" {
 #ifdef DOUBLE_PRECISION
 #  define kiss_fft_scalar celt_int32_t
 #  define kiss_twiddle_scalar celt_int32_t
+#  define KF_SUFFIX _celt_double
 #else
 #  define kiss_fft_scalar celt_int16_t
 #  define kiss_twiddle_scalar celt_int16_t
+#  define KF_SUFFIX _celt_single
 #endif
 #else
 # ifndef kiss_fft_scalar
 /*  default is float */
 #   define kiss_fft_scalar float
 #   define kiss_twiddle_scalar float
+#   define KF_SUFFIX _celt_single
 # endif
 #endif
 
 
 /* This adds a suffix to all the kiss_fft functions so we
    can easily link with more than one copy of the fft */
-#define KF_SUFFIX _celt_single
 #define CAT_SUFFIX(a,b) a ## b
 #define SUF(a,b) CAT_SUFFIX(a, b)
 
index 392f3eb..1f63365 100644 (file)
@@ -16,8 +16,10 @@ Redistribution and use in source and binary forms, with or without modification,
 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 COPYRIGHT OWNER 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"
+#ifndef SKIP_CONFIG_H
+#  ifdef HAVE_CONFIG_H
+#    include "config.h"
+#  endif
 #endif
 
 #include "os_support.h"
index c991c94..c8b3aa0 100644 (file)
 #include "config.h"
 #endif
 
+/*#include "_kiss_fft_guts.h"
+#include "kiss_fftr.h"*/
+#include "kfft_single.h"
+
 #include <stdio.h>
 #include <math.h>
 #include "pitch.h"
 #include "psy.h"
-#include "_kiss_fft_guts.h"
-#include "kiss_fftr.h"
+#include "os_support.h"
+
+kiss_fftr_cfg pitch_state_alloc(int max_lag)
+{
+   return kiss_fftr_alloc_celt_single(max_lag, 0, 0);
+}
+
+void pitch_state_free(kiss_fftr_cfg st)
+{
+   kiss_fft_free(st);
+}
+
+#ifdef FIXED_POINT
+static void normalise16(celt_word16_t *x, int len, celt_word16_t val)
+{
+   int i;
+   celt_word16_t maxval = 0;
+   for (i=0;i<len;i++)
+      maxval = MAX16(maxval, ABS16(x[i]));
+   if (maxval > val)
+   {
+      int shift = 0;
+      while (maxval > val)
+      {
+         maxval >>= 1;
+         shift++;
+      }
+      if (shift==0)
+         return;
+      for (i=0;i<len;i++)
+         x[i] = SHR16(x[i], shift);
+   } else {
+      int shift=0;
+      if (maxval == 0)
+         return;
+      val >>= 1;
+      while (maxval < val)
+      {
+         val >>= 1;
+         shift++;
+      }
+      if (shift==0)
+         return;
+      for (i=0;i<len;i++)
+         x[i] = SHL16(x[i], shift);
+   }
+}
+#else
+#define normalise16(x,len,val)
+#endif
+
+#define INPUT_SHIFT 15
 
 void find_spectral_pitch(kiss_fftr_cfg fft, struct PsyDecay *decay, const celt_sig_t *x, const celt_sig_t *y, const celt_word16_t *window, int overlap, int lag, int len, int C, int *pitch)
 {
    int c, i;
    celt_word32_t max_corr;
-   VARDECL(celt_word32_t *X);
-   VARDECL(celt_word32_t *Y);
+   VARDECL(celt_word16_t *X);
+   VARDECL(celt_word16_t *Y);
    VARDECL(celt_mask_t *curve);
    int n2;
    int L2;
@@ -59,7 +113,7 @@ void find_spectral_pitch(kiss_fftr_cfg fft, struct PsyDecay *decay, const celt_s
    SAVE_STACK;
    n2 = lag/2;
    L2 = len/2;
-   ALLOC(X, lag, celt_word32_t);
+   ALLOC(X, lag, celt_word16_t);
    ALLOC(curve, n2, celt_mask_t);
 
    bitrev = fft->substate->bitrev;
@@ -70,8 +124,8 @@ void find_spectral_pitch(kiss_fftr_cfg fft, struct PsyDecay *decay, const celt_s
    {
       for (i=0;i<L2;i++)
       {
-         X[2*bitrev[i]] += SHR32(x[C*(2*i)+c],1);
-         X[2*bitrev[i]+1] += SHR32(x[C*(2*i+1)+c],1);
+         X[2*bitrev[i]] += SHR32(x[C*(2*i)+c],INPUT_SHIFT);
+         X[2*bitrev[i]+1] += SHR32(x[C*(2*i+1)+c],INPUT_SHIFT);
       }
    }
    /* Applying the window in the bit-reverse domain. It's a bit weird, but it
@@ -83,6 +137,8 @@ void find_spectral_pitch(kiss_fftr_cfg fft, struct PsyDecay *decay, const celt_s
       X[2*bitrev[L2-i-1]] = MULT16_32_Q15(window[2*i+1], X[2*bitrev[L2-i-1]]);
       X[2*bitrev[L2-i-1]+1] = MULT16_32_Q15(window[2*i], X[2*bitrev[L2-i-1]+1]);
    }
+   normalise16(X, lag, 8192);
+   /*for (i=0;i<lag;i++) printf ("%d ", X[i]);printf ("\n");*/
    /* Forward real FFT (in-place) */
    kf_work((kiss_fft_cpx*)X, NULL, 1,1, fft->substate->factors,fft->substate, 1, 1, 1);
    kiss_fftr_twiddles(fft,X);
@@ -90,7 +146,7 @@ void find_spectral_pitch(kiss_fftr_cfg fft, struct PsyDecay *decay, const celt_s
    compute_masking(decay, X, curve, lag);
 
    /* Deferred allocation to reduce peak stack usage */
-   ALLOC(Y, lag, celt_word32_t);
+   ALLOC(Y, lag, celt_word16_t);
    for (i=0;i<lag;i++)
       Y[i] = 0;
    /* Sum all channels of the past audio and copy into Y in bit-reverse order */
@@ -98,10 +154,11 @@ void find_spectral_pitch(kiss_fftr_cfg fft, struct PsyDecay *decay, const celt_s
    {
       for (i=0;i<n2;i++)
       {
-         Y[2*bitrev[i]] += SHR32(y[C*(2*i)+c],1);
-         Y[2*bitrev[i]+1] += SHR32(y[C*(2*i+1)+c],1);
+         Y[2*bitrev[i]] += SHR32(y[C*(2*i)+c],INPUT_SHIFT);
+         Y[2*bitrev[i]+1] += SHR32(y[C*(2*i+1)+c],INPUT_SHIFT);
       }
    }
+   normalise16(Y, lag, 8192);
    /* Forward real FFT (in-place) */
    kf_work((kiss_fft_cpx*)Y, NULL, 1,1, fft->substate->factors,fft->substate, 1, 1, 1);
    kiss_fftr_twiddles(fft,Y);
@@ -113,7 +170,8 @@ void find_spectral_pitch(kiss_fftr_cfg fft, struct PsyDecay *decay, const celt_s
       celt_word32_t tmp;
       /*n = 1.f/(1e1+sqrt(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 = SIG_SCALING_1/sqrt(1+curve[i]);
+      /*printf ("%d %d ", X[2*i]*X[2*i]+X[2*i+1]*X[2*i+1], Y[2*i]*Y[2*i]+Y[2*i+1]*Y[2*i+1]);*/
+      n = 1.f/sqrt(1+curve[i]);
       /*printf ("%f ", n);*/
       /*n = 1.f/(1+curve[i]);*/
       tmp = X[2*i];
@@ -122,8 +180,11 @@ void find_spectral_pitch(kiss_fftr_cfg fft, struct PsyDecay *decay, const celt_s
    }
    /*printf ("\n");*/
    X[0] = X[1] = 0;
+   /*for (i=0;i<lag;i++) printf ("%d ", X[i]);printf ("\n");*/
+   normalise16(X, lag, 50);
    /* Inverse half-complex to real FFT gives us the correlation */
    kiss_fftri(fft, X, Y);
+   /*for (i=0;i<lag;i++) printf ("%d ", Y[i]);printf ("\n");*/
    /*for (i=0;i<C*lag;i++)
       printf ("%d %d\n", X[i], xx[i]);*/
    
index fe916b9..062f126 100644 (file)
@@ -41,6 +41,9 @@
 #include "kiss_fftr.h"
 #include "psy.h"
 
+kiss_fftr_cfg pitch_state_alloc(int max_lag);
+void pitch_state_free(kiss_fftr_cfg st);
+
 /** Find the optimal delay for the pitch prediction. Computation is
     done in the frequency domain, both to save time and to make it
     easier to apply psychoacoustic weighting */
index 294e2dc..ed8ec62 100644 (file)
@@ -126,7 +126,7 @@ static void spreading_func(struct PsyDecay *d, celt_word32_t *psd, celt_mask_t *
 }
 
 /* Compute a marking threshold from the spectrum X. */
-void compute_masking(struct PsyDecay *decay, celt_word32_t *X, celt_mask_t *mask, int len)
+void compute_masking(struct PsyDecay *decay, celt_word16_t *X, celt_mask_t *mask, int len)
 {
    int i;
    VARDECL(celt_word32_t *psd);
@@ -134,10 +134,10 @@ void compute_masking(struct PsyDecay *decay, celt_word32_t *X, celt_mask_t *mask
    SAVE_STACK;
    N=len/2;
    ALLOC(psd, N, celt_word32_t);
-   psd[0] = MULT16_16(ROUND(X[0],SIG_SHIFT), ROUND(X[0],SIG_SHIFT));
+   psd[0] = MULT16_16(X[0], X[0]);
    for (i=1;i<N;i++)
-      psd[i] = ADD32(MULT16_16(ROUND(X[i*2],SIG_SHIFT), ROUND(X[i*2],SIG_SHIFT)),
-                     MULT16_16(ROUND(X[i*2+1],SIG_SHIFT), ROUND(X[i*2+1],SIG_SHIFT)));
+      psd[i] = ADD32(MULT16_16(X[i*2], X[i*2]),
+                     MULT16_16(X[i*2+1], X[i*2+1]));
    /* TODO: Do tone masking */
    /* Noise masking */
    spreading_func(decay, psd, mask, N);
index 5180ce0..667bebd 100644 (file)
@@ -45,7 +45,7 @@ void psydecay_init(struct PsyDecay *decay, int len, celt_int32_t Fs);
 void psydecay_clear(struct PsyDecay *decay);
 
 /** Compute the masking curve for an input (DFT) spectrum X */
-void compute_masking(struct PsyDecay *decay, celt_word32_t *X, celt_mask_t *mask, int len);
+void compute_masking(struct PsyDecay *decay, celt_word16_t *X, celt_mask_t *mask, int len);
 
 /** Compute the masking curve for an input (MDCT) spectrum X */
 void compute_mdct_masking(struct PsyDecay *decay, celt_word32_t *X, celt_mask_t *mask, int len);