From 0058cbc913979291fcd57deb0b92d7a2e5a58fcf Mon Sep 17 00:00:00 2001 From: evilsocket Date: Tue, 26 May 2009 18:23:35 +0200 Subject: [PATCH] Experimental implementation of FDCT. --- Makefile | 16 +++- include/fdct.h | 53 +++++++++++ vocal.h => include/vocal.h | 0 dsp.c => src/dsp.c | 0 src/fdct.c | 181 +++++++++++++++++++++++++++++++++++++ main.c => src/main.c | 83 ++--------------- utils.c => src/utils.c | 0 7 files changed, 257 insertions(+), 76 deletions(-) create mode 100644 include/fdct.h rename vocal.h => include/vocal.h (100%) rename dsp.c => src/dsp.c (100%) create mode 100644 src/fdct.c rename main.c => src/main.c (69%) rename utils.c => src/utils.c (100%) diff --git a/Makefile b/Makefile index 28c6153..16b813a 100644 --- a/Makefile +++ b/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 diff --git a/include/fdct.h b/include/fdct.h new file mode 100644 index 0000000..f72e4a6 --- /dev/null +++ b/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 +#include +#include + +#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 ); diff --git a/vocal.h b/include/vocal.h similarity index 100% rename from vocal.h rename to include/vocal.h diff --git a/dsp.c b/src/dsp.c similarity index 100% rename from dsp.c rename to src/dsp.c diff --git a/src/fdct.c b/src/fdct.c new file mode 100644 index 0000000..e9446fe --- /dev/null +++ b/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); +} diff --git a/main.c b/src/main.c similarity index 69% rename from main.c rename to src/main.c index 9d42a3b..abc2901 100644 --- a/main.c +++ b/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,8 +64,7 @@ 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); + memset(buf, 0x80, TOTSIZE ); if (!fname) { fname = (char*) malloc(0x100); @@ -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); - - 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); - } + fdct_init( &dct, buf ); - free(neutral); - free(buf); + sum = fdct_sum( &dct, 36927.96, NULL/*neutral*/, TOTSIZE ); + + fdct_release(&dct); t2 = time((unsigned) NULL); + printf ("DCT computing: done in %u seconds\n", (unsigned int) (t2-t1)); - - for (i=0; i