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.

#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;
 
  }
 
 
}
 
Except where otherwise noted, content on this wiki is licensed under the following license: CC Attribution-Noncommercial-Share Alike 4.0 International
Recent changes RSS feed Donate Powered by PHP Valid XHTML 1.0 Valid CSS Driven by DokuWiki