Hanging on IRC, we decided to try how fast can I write a DTMF decoder with library FFT. I estimated 30 minutes. Finally, it took me exactly one hour, and the code quality is terrible.
Yes, I know that computing these 8 bins from definition would be faster and would not require fftw.
Now the whole world sees how I code for competitions where only right output, not the code quality, is evaluated.
#include <stdlib.h> #include <stdio.h> #include <inttypes.h> #include <string.h> #include <unistd.h> #include <math.h> #include <fftw3.h> #include <fcntl.h> #include <sys/mman.h> #include <sys/types.h> #include <unistd.h> #include <sys/stat.h> #include <sys/ioctl.h> #include <string.h> #include <linux/fs.h> double hamming(int i, int length) { double a, b, w, N1; a = 25.0/46.0; b = 21.0/46.0; N1 = (double)(length-1); w = a - b*cos(2*i*M_PI/N1); return w; } int main(int argc, char **argv) { int i; float* acc; float* fftbuf; int fftsize = 512; acc = (float*) calloc(fftsize * 2, sizeof(float)); fftbuf = (float*) calloc(fftsize * 2, sizeof(float)); /* Setup FFTW */ int N = fftsize; fftwf_complex *fftw_in, *fftw_out; fftwf_plan p; fftw_in = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * N); fftw_out = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * N); p = fftwf_plan_dft_1d(N, fftw_in, fftw_out, FFTW_FORWARD, FFTW_ESTIMATE); float * win = malloc(fftsize * sizeof(float)); for(i=0; i<fftsize; i++) { win[i] = hamming(i, fftsize); } int16_t * data; int fd = open("/tmp/nahravka.wav", O_RDONLY); size_t dsize = 168204; data = (int16_t*)mmap(NULL, dsize, PROT_READ, MAP_PRIVATE, fd, 0); int samples = dsize/2; /* compute smooth func */ int16_t * mx = malloc(sizeof(int16_t) * dsize); for(int i = 0; i<samples-50; i++) { int max = 0; for(int j = i; j<i+30; j++) { if(abs(data[j]) > max) { max = abs(data[j]); } } mx[i] = max; /*if(i%20 == 1) { printf("max %i\n", max); }*/ } #define sigth1 15000 #define sigth2 10000 #define sigstart 32000 //printf("HERE\n"); int pos = sigstart; while(1) { /* find start of burst */ while(pos++) { // printf("pos %i mx %i\n", pos, mx[pos]); if(mx[pos] > sigth1) { break; } } /* find end of burst */ int endpos = pos+10; while(endpos++) { // printf("endpos %i mx %i\n", endpos, mx[endpos]); if(mx[endpos] < sigth2) { break; } } /* align */ int center = (pos+endpos)/2; int dur = endpos-pos; int start = center-fftsize/2; //printf("len = %i (%i - %i), start %i\n", endpos-pos, pos, endpos, start); /* Read data */ for(i=0; i<fftsize; i++) { float sample = (float) ( data[i+start] ); fftbuf[i*2] = sample * win[i]; fftbuf[i*2+1] = 0.0; } // Transform memcpy(fftw_in, fftbuf, 2*fftsize * sizeof(float)); fftwf_execute(p); float* fout = (float*)fftw_out; // Read output float* rout = (float*)malloc(sizeof(float)*fftsize); for(int i = 0; i<fftsize; i++) { float orig = hypot(fftbuf[i*2],fftbuf[i*2+1]); float ret = hypot(fout[i*2], fout[i*2+1]); rout[i] = ret; //printf("%f %f\n", orig, ret); } /* find two maximums */ int max1 = 8; int max2 = 8; rout[511] = 0.0; float dtmf[] = { 1209, 1336, 1477, 1633, 697, 770, 852, 941 }; int bins[9]; for(int i = 0; i<8; i++) { bins[i] = (float)fftsize * dtmf[i]/8000.0; //printf("bin %i = %f\n", bins[i], rout[bins[i]]); } bins[8] = 511; /* horiz */ for(int i = 0; i<4; i++) { if(rout[bins[i]] > rout[bins[max1]]) { max1 = i; } } /* vert */ for(int i = 4; i<8; i++) { //printf("%i cmp %f %f\n", i, rout[bins[i]], rout[bins[max2]]); if(rout[bins[i]] > rout[bins[max2]]) { max2 = i; } } max2 -= 4; //printf("%i %i\n", max1, max2); int res[4][4] = {{1, 2, 3, 'A'}, {4, 5, 6, 'B'}, {7, 8, 9, 'C'}, {'*', 0, '#', 'D'} }; int rt = res[max2][max1]; if(rt>10) { printf(" *** %c ***\n", rt); } else { printf(" *** %i ***\n", rt); } pos = endpos+10; } }