Experimental implementation of FDCT.

Makefile

@@ -1,9 +1,19 @@
+CC = gcc
+SRCDIR = src/
+INCDIR = include/
+CFLAGS = -I$(INCDIR) -O3 -funroll-loops -w -ffast-math -fno-stack-protector -ffunction-sections -funsafe-math-optimizations -fno-trapping-math
+LFLAGS = -lm
+
 all:
-	gcc -Wall -pedantic -o vocal main.c dsp.c utils.c -lm -O3 -funroll-loops -w -ffast-math -fno-stack-protector -ffunction-sections -funsafe-math-optimizations -fno-trapping-math
+	$(CC) -c $(CFLAGS) $(SRCDIR)utils.c
+	$(CC) -c $(CFLAGS) $(SRCDIR)dsp.c
+	$(CC) -c $(CFLAGS) $(SRCDIR)fdct.c
+	
+	$(CC) $(CFLAGS) $(LFLAGS) $(SRCDIR)main.c *.o -o voxifera 
 
 clean:
-	rm vocal
+	rm -f voxifera *.o
 
 install:
 	mkdir -p /usr/local/bin
-	cp vocal /usr/local/bin
+	cp voxifera /usr/local/bin

include/fdct.h

@@ -0,0 +1,53 @@
+/***************************************************************************
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ *   This program 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 General Public License for more details.                          *
+ *                                                                         *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+
+// Very fast DCT (discrete cosine transform) implementation based on IEEE signal proc, 1992 #8 .
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+ 
+#define FDCT_TOTSIZE 16384 /* !!! QUESTA È LA GRANDEZZA DEL BUFFER --DEVE-- ESSERE UNA POTENZA DI 2 !!! */
+#define FDCT_POW 14     /* ln( FDCT_TOTSIZE ) */
+#define FDCT_TSO2 8192 /* FDCT_TOTSIZE / 2 */
+#define FDCT_TSO4 4096     /* FDCT_TOTSIZE / 4 */
+#define FDCT_2OTS 0.00012207 /* 2 / FDCT_TOTSIZE */
+#define FDCT_2OTS_SQ 0.011048543 /* sqrt( 2 / FDCT_TOTSIZE ) */
+#define FDCT_FACTOR 0.7071067814
+
+// fast DCT main structure 
+typedef struct {
+  double * transform;
+  double * ctable;
+}
+fdct_t;
+
+// initialize the structure with data
+void fdct_init( fdct_t *dct, unsigned char *data );
+// reverse data 'til size-th byte
+inline void fdct_bitrev( double *data, unsigned int size );
+// complete data reversing
+inline void fdct_scramble( double *data );
+// compute cosine coefficents
+inline void fdct_butterflies( fdct_t *dct );
+// compute FDCT sums on data
+inline void fdct_sums(double *data);
+// compute FDCT
+inline void fdct_compute( fdct_t *dct );
+// compute FDCT vector sum, where k is the normalizing factor and neutral an optional mask to subtract
+inline double fdct_sum( fdct_t *dct, double div, unsigned char *neutral, unsigned int neutral_size  );
+// release FDCT structure
+void fdct_release( fdct_t *dct );

src/fdct.c

@@ -0,0 +1,181 @@
+/***************************************************************************
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ *   This program 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 General Public License for more details.                          *
+ *                                                                         *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+#include "fdct.h"
+
+void fdct_init( fdct_t *dct, unsigned char *data ){
+  int i, group, base, item, nitems;
+    double factor;
+  
+  /* alloc transform vector and fill with data */
+  dct->transform = (double *)calloc( FDCT_TOTSIZE, sizeof(double) );
+  for( i = 0; i < FDCT_TOTSIZE; i++ ){
+    dct->transform[i] = (double)data[i];
+  }
+  
+  /* alloc and fill cosines table */
+  dct->ctable = (double *)calloc( FDCT_TOTSIZE, sizeof(double) );  
+  for( i = 0; i <= FDCT_TSO2 - 1; i++ ){
+    dct->ctable[FDCT_TSO2 + i] = 4 * i + 1;    
+  }
+  
+  for( group = 1; group <= FDCT_POW - 1; group++ ){
+    nitems = base = 1 << (group - 1);
+    factor = 1.0 * (1 << (FDCT_POW - group));
+    for( item = 1; item <= nitems; item++ ){
+      dct->ctable[base + item - 1] = factor * dct->ctable[FDCT_TSO2 + item - 1];
+    }
+  }
+ 
+  for( i = 1; i <= FDCT_TOTSIZE - 1; i++ ){
+    dct->ctable[i] = 1.0 / (2.0 * cos( dct->ctable[i] * M_PI / (2.0 * FDCT_TOTSIZE) ));
+  }
+  
+  return dct;
+}
+ 
+inline void fdct_bitrev( double *data, unsigned int size ){
+  int i, m,
+    j = 1,
+    hsize = size >> 1;
+  double tmp;
+ 
+  if( size <= 2 ){ return; }
+ 
+  for( i = 1; i <= size; i++ ){
+    if(i < j){
+      tmp = data[j - 1];
+      data[j - 1] = data[i - 1];
+      data[i - 1] = tmp;
+    }
+    m = hsize;
+    while(j > m){
+      j = j - m;
+      m = (m + 1) >> 1;
+    }
+    j = j + m;
+  }
+}
+ 
+inline void fdct_scramble( double *data ){
+  double tmp;
+  int i, ii1, ii2;
+ 
+  fdct_bitrev( data, FDCT_TOTSIZE );
+  fdct_bitrev( &data[0], FDCT_TSO2 );
+  fdct_bitrev( &data[FDCT_TSO2], FDCT_TSO2);
+ 
+  ii1 = FDCT_TOTSIZE - 1;
+  ii2 = FDCT_TSO2;
+  for( i = 0; i <= FDCT_TSO4 - 1; i++, ii1--, ii2++ ){
+    tmp = data[ii1];
+    data[ii1] = data[ii2];
+    data[ii2] = tmp;
+  }
+}
+ 
+inline void fdct_butterflies( fdct_t *dct ){
+  int stage, ii1, ii2,
+    butterfly,
+    ngroups,
+    group,
+    wingspan,
+    inc,
+    ptr;
+  double cosine, value;
+ 
+  for( stage = FDCT_POW; stage >= 1; stage-- ){
+    ngroups = 1 << (FDCT_POW - stage);
+    wingspan = 1 << (stage - 1);
+    inc = wingspan << 1;
+        
+    for( butterfly = 1; butterfly <= wingspan; butterfly++ ){
+      cosine = dct->ctable[ wingspan + butterfly - 1 ];
+      ptr = 0;
+      
+      for( group = 1; group <= ngroups; group++ ){
+        ii1 = ptr + butterfly - 1;
+        ii2 = ii1 + wingspan;
+        value = dct->transform[ii2];
+        dct->transform[ii2] = cosine *(dct->transform[ii1] - value);
+        dct->transform[ii1] = dct->transform[ii1] + value;
+        ptr += inc;
+      }
+    }
+  }
+}
+ 
+inline void fdct_sums(double *data){
+  int stepsize,
+    stage, i,
+    nthreads,
+    thread,
+    step,
+    nsteps;
+ 
+  for( stage = FDCT_POW - 1; stage >= 1; stage-- ){
+    nthreads = 1 << (stage - 1);
+    stepsize = nthreads << 1;
+    nsteps = (1 << (FDCT_POW - stage)) - 1;
+    
+    for( thread = 1; thread <= nthreads; thread++ ){
+      i = nthreads + thread - 1;
+      for( step = 1; step <= nsteps; step++ ){
+        data[i] += data[i + stepsize];
+        i += stepsize;
+      }
+    }
+  }
+}
+ 
+inline void fdct_compute( fdct_t *dct ){
+  int i;
+  
+  fdct_scramble( dct->transform );
+    fdct_butterflies( dct );
+    fdct_bitrev( dct->transform, FDCT_TOTSIZE );
+    fdct_sums(dct->transform);
+    
+  dct->transform[0] *= FDCT_FACTOR;
+  for( i = 0; i <= FDCT_TOTSIZE - 1; i++ ){
+    dct->transform[i] *= FDCT_2OTS_SQ;
+  }    
+}
+ 
+inline double fdct_sum( fdct_t *dct, double div, unsigned char *neutral, unsigned int neutral_size ){
+  int i;
+  double sum = 0.0;
+  
+  fdct_compute( dct );
+  
+  if( neutral ){
+    for( i = 0; i < FDCT_TOTSIZE; i++ ){
+      sum += dct->transform[i] - neutral[i % neutral_size];
+    }
+  }
+  else{
+    for( i = 0; i < FDCT_TOTSIZE; i++ ){
+      sum += dct->transform[i];
+    }
+  }
+  
+    return fabs(sum / div);
+}
+ 
+void fdct_release( fdct_t *dct ){
+  free(dct->ctable);
+  free(dct->transform);
+} 

src/main.c

@@ -16,43 +16,10 @@
  ***************************************************************************/
 
 #include "vocal.h"
+#include "fdct.h"
 
 typedef enum  { capture, append } mode;
 
-/* tabella dei fattoriali di 1/n!, dove 0 <= n < 30 */
-static double ftable[]  = { 0.0, 1.0, 0.5, 0.166667, 0.0416667, 0.00833333, 0.00138889, 0.000198413, 
-						    2.48016e-05, 2.75573e-06, 2.75573e-07, 2.50521e-08, 2.08768e-09, 5.17584e-10, 
-						    7.81894e-10, 4.98925e-10, 4.98955e-10, -3.46594e-09, -1.11305e-09, 9.12062e-09, 
-						   -4.75707e-10, -8.3674e-10, -1.91309e-09, 1.15948e-09, -1.28875e-09, 4.81654e-10, 
-						   -5.39409e-10, 6.73499e-10, -7.26886e-10, -8.05468e-10 };
-/* intervallo in cui posizionare x */
-static double offset[]  = { 0.0, M_PI };
-/* sviluppo in serie di Taylor della funzione coseno 
- * 
- * x     : valore per cui calcolare il coseno .
- * terms : numero di termini da usare nello sviluppo dove 0 <= terms < 30 e terms è un numero pari .
- */
-inline double taylor_cosine( double x, int terms ) {
-	
-    int i        = terms,
-	    quadrant = x * TWO_O_M_PI; /* 0, 1, 2 o 3 */
-	double x2, r;
-	
-	/* tariamo x in modo tale che −π < x < π */
-    x         = x - quadrant * H_M_PI;
-    quadrant += 1;
-    x         = offset[ (quadrant >> 1) & 1 ] - x;
-    x2        = -(x*x);
-	/* eseguo lo sviluppo in serie */
-    r = ftable[i] * x2;
-    for( i -= 2; i > 0; i -= 2 ){
-    	r += ftable[i];
-        r *= x2;
-    }
-
-    return r + 1;
-}
-
 void helpandusage()  {
 	fprintf (stderr,
 			"---=== BlackLight's & evilsocket's vocal recognition v.0.1b ===---\n"
@@ -72,10 +39,9 @@ int main (int argc, char **argv)  {
 	char *line, *tmp, *cmd = NULL, *fname = NULL, *device = NULL;
 	char **match = NULL;
 	
-	double deviance;
-	double *neutral, *trans;
-	double value, sum = 0.0;
-	double t, v;
+	fdct_t dct;
+	double sum, deviance, value;
+	double *neutral;
 	
 	FILE *fp;
 	time_t t1, t2;
@@ -98,7 +64,6 @@ int main (int argc, char **argv)  {
 
 	buf     = (u8*) malloc(TOTSIZE);
 	neutral = (double*) malloc(TOTSIZE*sizeof(double));
-	trans   = (double*) malloc(TOTSIZE*sizeof(double));
 	memset(buf, 0x80, TOTSIZE );
 
 	if (!fname)  {
@@ -131,43 +96,15 @@ int main (int argc, char **argv)  {
 	t1 = time((unsigned) NULL);
 	
 	
-	/* porto fuori il loop per i = 0 così da evitare controlli su t (che dipende da i) 
-	 * e j nel secondo loop più consistente . */
-	v = 0;
-	for( j = 0; j < TOTSIZE; ++j ){
-		v += buf[j];
-	}
-	v *= M_COEFF;
-	v -= neutral[0];
-	v = (v >= 0) ? v : -v;
-	trans[0] = log(v);
+	fdct_init( &dct, buf );
 
-	double init = buf[0];
-	for( i = 1; i < TOTSIZE; ++i ) {
-		t = M_FACTOR * i;
-		v = init;
-		/* dato che ho inizializzato v a buf[0] e sono sicuro che t != 0, 
-		 * posso partire da j = 1 . */
-		for( j = 1; j < TOTSIZE; ++j ){
-			v += buf[j] * taylor_cosine( t * j, 10 );
-		}
-		v *= M_COEFF;
-		v -= neutral[i];
-		v = (v >= 0) ? v : -v;
-		trans[i] = log(v);
-	}
+	sum = fdct_sum( &dct, 36927.96, NULL/*neutral*/, TOTSIZE );
 
-	free(neutral);
-	free(buf);
+	fdct_release(&dct);
 
 	t2 = time((unsigned) NULL);
+	
 	printf ("DCT computing: done in %u seconds\n", (unsigned int) (t2-t1));
-
-	for (i=0; i<TOTSIZE; i++)
-		sum += trans[i];
-	sum = exp(sum/100000.0);
-	free(trans);
-
 	printf ("vocal sequence successfully acquired, Fourier sum = %g\n", sum);
 
 	if (m == append)  {