====== FFTW demo - DTMF decoder ======

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. 
<code c>
#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;

  }


}

</code>