Added pitch analysis. Doesn't crash, but otherwise untested.

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"

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

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

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;
 }
 

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);
 
 
 

libcelt/fftwrap.c

+/* 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

libcelt/fftwrap.h

+/* 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

libcelt/pitch.c

+/**
+   @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);
+}

libcelt/pitch.h

+/**
+   @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

libcelt/smallft.h

+/********************************************************************
+ *                                                                  *
+ * 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

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;
 }