use replacement if complex.h not found
authorPaul Brossier <piem@altern.org>
Fri, 29 Sep 2006 14:19:55 +0000 (14:19 +0000)
committerPaul Brossier <piem@altern.org>
Fri, 29 Sep 2006 14:19:55 +0000 (14:19 +0000)
use replacement if complex.h not found

configure.ac
src/fft.c
src/onsetdetection.c

index e0cfdcce6d859fb9eb0f13dbb3753739108a1840..ee4786251507252c4ece39680fc995f58c057e05 100644 (file)
@@ -105,8 +105,14 @@ AC_CHECK_LIB(pthread, pthread_create)
 dnl Check for header files
 AC_HEADER_STDC
 AC_CHECK_HEADERS([string.h stdlib.h stdio.h math.h errno.h stdarg.h unistd.h signal.h],,)
-AC_CHECK_HEADERS(complex.h,,AC_MSG_ERROR([Ouch! missing complex.h header]))
 AC_CHECK_HEADERS(fftw3.h  ,,AC_MSG_ERROR([Ouch! missing fftw3.h header]))
+AC_ARG_ENABLE(complex,
+  AC_HELP_STRING([--enable-complex],[compile with complex.h [[default=auto]]]),
+  [with_complex=$enableval],
+  with_complex="yes")
+if test "$with_complex" = "yes"; then
+  AC_CHECK_HEADERS(complex.h,,AC_MSG_WARN([Ouch! missing complex.h header]))
+fi
 
 dnl Check for __VAR_ARGS__ support
 AC_CACHE_CHECK(for C99 __VA_ARGS__ macro,
index 01e3783053107310a39d7b60d72c79bf91cc2361..37c7cb900a6bebe0e29e891636ebd7a3e7d458df 100644 (file)
--- a/src/fft.c
+++ b/src/fft.c
@@ -28,6 +28,7 @@
 #define fftw_execute           fftwf_execute
 #define fftw_plan_dft_r2c_1d   fftwf_plan_dft_r2c_1d
 #define fftw_plan_dft_c2r_1d   fftwf_plan_dft_c2r_1d
+#define fftw_plan_r2r_1d      fftwf_plan_r2r_1d
 #define fftw_plan              fftwf_plan
 #define fftw_destroy_plan      fftwf_destroy_plan
 #endif
@@ -46,6 +47,8 @@ struct _aubio_fft_t {
        fftw_plan       pfw, pbw;
 };
 
+static void aubio_fft_getspectrum(fft_data_t * spectrum, smpl_t *norm, smpl_t * phas, uint_t size);
+
 aubio_fft_t * new_aubio_fft(uint_t size) {
        aubio_fft_t * s = AUBIO_NEW(aubio_fft_t);
        /* allocate memory */
@@ -53,8 +56,13 @@ aubio_fft_t * new_aubio_fft(uint_t size) {
        s->out      = AUBIO_ARRAY(real_t,size);
        s->specdata = (fft_data_t*)fftw_malloc(sizeof(fft_data_t)*size);
        /* create plans */
+#ifdef HAVE_COMPLEX_H
        s->pfw = fftw_plan_dft_r2c_1d(size, s->in,  s->specdata, FFTW_ESTIMATE);
        s->pbw = fftw_plan_dft_c2r_1d(size, s->specdata, s->out, FFTW_ESTIMATE);
+#else
+       s->pfw = fftw_plan_r2r_1d(size, s->in,  s->specdata, FFTW_R2HC, FFTW_ESTIMATE);
+       s->pbw = fftw_plan_r2r_1d(size, s->specdata, s->out, FFTW_HC2R, FFTW_ESTIMATE);
+#endif
        return s;
 }
 
@@ -88,17 +96,55 @@ void aubio_fft_rdo(const aubio_fft_t * s,
        for (i=0;i<size;i++) data[i] = s->out[i]*renorm;
 }
 
+#ifdef HAVE_COMPLEX_H
 
 void aubio_fft_getnorm(smpl_t * norm, fft_data_t * spectrum, uint_t size) {
        uint_t i;
-       for (i=0;i<size;i++) norm[i] = ABSC(spectrum[i]);
+       for (i=0;i<size/2+1;i++) norm[i] = ABSC(spectrum[i]);
+       //for (i=0;i<size/2+1;i++) AUBIO_DBG("%f\n", norm[i]);
 }
 
 void aubio_fft_getphas(smpl_t * phas, fft_data_t * spectrum, uint_t size) {
        uint_t i;
-       for (i=0;i<size;i++) phas[i] = ARGC(spectrum[i]);
+       for (i=0;i<size/2+1;i++) phas[i] = ARGC(spectrum[i]);
+       //for (i=0;i<size/2+1;i++) AUBIO_DBG("%f\n", phas[i]);
+}
+
+void aubio_fft_getspectrum(fft_data_t * spectrum, smpl_t *norm, smpl_t * phas, uint_t size) {
+  uint_t j;
+  for (j=0; j<size/2+1; j++) {
+    spectrum[j]  = CEXPC(I*phas[j]);
+    spectrum[j] *= norm[j];
+  }
 }
 
+#else
+
+void aubio_fft_getnorm(smpl_t * norm, fft_data_t * spectrum, uint_t size) {
+       uint_t i;
+  norm[0] = -spectrum[0];
+       for (i=1;i<size/2+1;i++) norm[i] = SQRT(SQR(spectrum[i]) + SQR(spectrum[size-i]));
+       //for (i=0;i<size/2+1;i++) AUBIO_DBG("%f\n", norm[i]);
+}
+
+void aubio_fft_getphas(smpl_t * phas, fft_data_t * spectrum, uint_t size) {
+       uint_t i;
+  phas[0] = PI;
+       for (i=1;i<size/2+1;i++) phas[i] = atan2f(spectrum[size-i] , spectrum[i]);
+       //for (i=0;i<size/2+1;i++) AUBIO_DBG("%f\n", phas[i]);
+}
+
+void aubio_fft_getspectrum(fft_data_t * spectrum, smpl_t *norm, smpl_t * phas, uint_t size) {
+  uint_t j;
+  for (j=0; j<size/2+1; j++) {
+    spectrum[j]       = norm[j]*COS(phas[j]);
+  }
+  for (j=1; j<size/2+1; j++) {
+    spectrum[size-j]  = norm[j]*SIN(phas[j]);
+  }
+}
+
+#endif
 
 /* new interface aubio_mfft */
 struct _aubio_mfft_t {
@@ -127,19 +173,16 @@ void aubio_mfft_do (aubio_mfft_t * fft,fvec_t * in,cvec_t * fftgrain){
         for (i=0; i < fft->channels; i++) {
                 aubio_fft_do (fft->fft,in->data[i],fft->spec[i],fft->winsize);
                 /* put norm and phase into fftgrain */
-                aubio_fft_getnorm(fftgrain->norm[i], fft->spec[i], fft->winsize/2+1);
-                aubio_fft_getphas(fftgrain->phas[i], fft->spec[i], fft->winsize/2+1);
+                aubio_fft_getnorm(fftgrain->norm[i], fft->spec[i], fft->winsize);
+                aubio_fft_getphas(fftgrain->phas[i], fft->spec[i], fft->winsize);
         }
 }
 
 /* execute inverse fourier transform */
 void aubio_mfft_rdo(aubio_mfft_t * fft,cvec_t * fftgrain, fvec_t * out){
-        uint_t i=0,j;
+        uint_t i=0;
         for (i=0; i < fft->channels; i++) {
-                for (j=0; j<fft->winsize/2+1; j++) {
-                        fft->spec[i][j]  = CEXPC(I*aubio_unwrap2pi(fftgrain->phas[i][j]));
-                        fft->spec[i][j] *= fftgrain->norm[i][j];
-                }
+                aubio_fft_getspectrum(fft->spec[i],fftgrain->norm[i],fftgrain->phas[i],fft->winsize);
                 aubio_fft_rdo(fft->fft,fft->spec[i],out->data[i],fft->winsize);
         }
 }
index 71f3c921c757a17e6798984f8461bf1560769979..2514c640bb6de1c77a84c64e547b7558ed128e88 100644 (file)
@@ -76,12 +76,22 @@ void aubio_onsetdetection_complex (aubio_onsetdetection_t *o, cvec_t * fftgrain,
                                        fftgrain->phas[i][j]
                                        -2.0*o->theta1->data[i][j]+
                                        o->theta2->data[i][j]);
+#ifdef HAVE_COMPLEX_H
                        o->meas[j] = fftgrain->norm[i][j]*CEXPC(I*o->dev1->data[i][j]);
                        /* sum on all bins */
                        onset->data[i][0]                += //(fftgrain->norm[i][j]);
                                        SQRT(SQR( REAL(o->oldmag->data[i][j]-o->meas[j]) )
                                                +  SQR( IMAG(o->oldmag->data[i][j]-o->meas[j]) )
                                                );
+#else
+                       o->meas[j]             = (fftgrain->norm[i][j])*COS(o->dev1->data[i][j]);
+                       o->meas[(nbins-1)*2-j] = (fftgrain->norm[i][j])*SIN(o->dev1->data[i][j]);
+                       /* sum on all bins */
+                       onset->data[i][0]                += //(fftgrain->norm[i][j]);
+                                       SQRT(SQR( (o->oldmag->data[i][j]-o->meas[j]) )
+                                               +  SQR( (-o->meas[(nbins-1)*2-j]) )
+                                               );
+#endif
                        /* swap old phase data (need to remember 2 frames behind)*/
                        o->theta2->data[i][j] = o->theta1->data[i][j];
                        o->theta1->data[i][j] = fftgrain->phas[i][j];
@@ -199,6 +209,7 @@ new_aubio_onsetdetection (aubio_onsetdetection_type type,
                uint_t size, uint_t channels){
        aubio_onsetdetection_t * o = AUBIO_NEW(aubio_onsetdetection_t);
        uint_t rsize = size/2+1;
+  uint_t i;
        switch(type) {
                /* for both energy and hfc, only fftgrain->norm is required */
                case aubio_onset_energy: 
@@ -209,7 +220,8 @@ new_aubio_onsetdetection (aubio_onsetdetection_type type,
                case aubio_onset_complex:
                        o->oldmag = new_fvec(rsize,channels);
                        /** bug: must be complex array */
-                       o->meas = AUBIO_ARRAY(fft_data_t,size);
+                       o->meas = AUBIO_ARRAY(fft_data_t,size+1);
+                       for (i=0; i<size+1; i++) o->meas[i] = 0;
                        o->dev1  = new_fvec(rsize,channels);
                        o->theta1 = new_fvec(rsize,channels);
                        o->theta2 = new_fvec(rsize,channels);