/* phoneman -- telephone manager
   A.J. Fisher	 January 1996
   DTMF recognition routines */

#include <stdio.h>	/* ??? */
#include <fishaudio.h>

#include "complex.h"
#include "phoneman.h"

#define TIMEOUT	    10			/* secs */
#define FFTDEPTH    8
#define FFTLENGTH   (1 << FFTDEPTH)
#define ALPHA	    0.54		/* for Hamming window */
#define THRESHOLD   1e14
#define NUMFREQS    9

#define DIG_INVALID (-1)
#define DIG_NOTONE  (-2)
#define DIG_TIMEOUT (-3)

static float freqs[NUMFREQS] =
  {  697.0,  770.0,  852.0,  941.0,	/* row freqs	*/
    1209.0, 1336.0, 1477.0, 1633.0,	/* column freqs */
    1100.0,				/* fax CNG tone */
  };

static uchar shufftab[FFTLENGTH];
static complex circle[FFTLENGTH/2 + 1]; /* extra slot needed for window computation */
static schar digits[1 << NUMFREQS];
static tone *dt1, *dt2;
static int time;

static uchar rotate(uchar, int);
static int getdigit(bool);
static float window(int);
static float powerat(complex*, float);
static void debugdigit(int);
static void dofft(complex*);
static void shuffle(complex*, complex*);
static void sub_fft(complex*, complex*, complex*, int);

inline float power(complex z) { return (z.re*z.re) + (z.im*z.im); }


global void initdtmf()	/* called at startup to initialize tables */
  { dt1 = new tone(350.0); dt2 = new tone(450.0);   /* 2 components of dial tone */
    for (int i=0; i < FFTLENGTH; i++) shufftab[i] = rotate(i, FFTDEPTH);
    for (int i=0; i <= FFTLENGTH/2; i++)
      { double w = TWOPI / FFTLENGTH;
	circle[i] = complex(cos(i*w), -sin(i*w));
      }
    digits[0] = DIG_NOTONE;
    for (int i=1; i < (1 << NUMFREQS); i++) digits[i] = DIG_INVALID;
    for (int i=0; i < 17; i++)
      { static ushort dt[17] =
	  { 0x028, 0x011, 0x021, 0x041, 0x012, 0x022, 0x042, 0x014,	/* 01234567 */
	    0x024, 0x044, 0x018, 0x048, 0x081, 0x082, 0x083, 0x084,	/* 89*#ABCD */
	    0x100,							/* fax CNG  */
	  };
	digits[dt[i]] = i;
      }
  }

static uchar rotate(uchar n, int d)
  { if (d > 1)
      { int odd = n & 1;
	n = rotate(n >> 1, d-1);
	if (odd) n |= (1 << (d-1));
      }
    return n;
  }

global void tidydtmf()
  { delete dt1; /* just to be tidy */
    delete dt2;
  }

global int getdtmf(bool first)
  { time = 0;
    putc('[', stderr);		/* ??? */
    int dig1, dig2;
    dig2 = getdigit(first);
    while (dig2 >= 0)
      { /* wait for silence or timeout */
	dig1 = dig2;
	dig2 = getdigit(first);
      }
    until ((dig2 >= 0 && dig2 == dig1) || dig2 == DIG_TIMEOUT)
      { /* wait for 2 identical tones, or timeout */
	dig1 = dig2;
	dig2 = getdigit(first);
      }
    fprintf(stderr, "]\n");     /* ??? */
    return dig2;
  }

static int getdigit(bool first)
  { if (first)
      { /* put dial tone in output buffer */
	while ((Audio -> ocount() + Audio -> icount()) < 3*TONELEN)
	  { for (int i=0; i < TONELEN; i++)
	      { int x = mu_expand(dt1 -> vec[i]) + mu_expand(dt2 -> vec[i]);
		Audio -> write(x << 7); /* 200 ms */
	      }
	  }
      }
    complex vec[FFTLENGTH];
    for (int j=0; j < FFTLENGTH; j++)
      { vec[j] = Audio -> read() * window(j);
	if (time++ > TIMEOUT*SAMPLERATE) return DIG_TIMEOUT;
      }
    dofft(vec);
    ushort result = 0;
    for (int j = NUMFREQS-1; j >= 0; j--)
      { float f = freqs[j];
	float p = powerat(vec, f) - powerat(vec, 2.0*f);  /* fundamental minus 2nd harmonic */
	result = (result << 1) | (p >= THRESHOLD);
      }
    unless (first) result &= ~0x100;	/* recognize fax CNG tone on first digit only */
    int d = digits[result];		/* return 0-15, or 16 if fax tone */
    debugdigit(d);
    return d;
  }

static float window(int j)
  { /* Hamming window */
    float win = (j < FFTLENGTH/2) ? circle[FFTLENGTH/2 - j].re : circle[j - FFTLENGTH/2].re;
    return ALPHA + (1.0-ALPHA) * win;
  }

static float powerat(complex *vec, float f)
  { float bn = f * (float) FFTLENGTH / (float) SAMPLERATE;  /* bin number */
    int n = (int) bn; float r = bn - (float) n;
    return (1.0-r) * power(vec[n]) + r * power(vec[n+1]);
  }

static void debugdigit(int d)
  { if (d < -2 || d > 16) giveup("Bug: ilgl digit %d", d);
    static char *debugstr = ".?0123456789*#ABCDf";
    putc(debugstr[d+2], stderr);
  }

/* fft routine generated by: /usr/fisher/mipsbin/mkfft -d 8 -g dofft
   ... and then massaged */

static void dofft(complex *data)
  { complex temp1[FFTLENGTH], temp2[FFTLENGTH];
    shuffle(data, temp1);
    sub_fft(temp1, data, temp2, FFTDEPTH);
  }

static void shuffle(complex *in, complex *out)
  { for (int i=0; i < FFTLENGTH; i++) out[shufftab[i]] = in[i];
  }

static void sub_fft(complex *in, complex *out, complex *temp, int d)
  { if (d > 0)
      { int h = 1 << (d-1), t = 1 << (FFTDEPTH-d);
	sub_fft(&in[0], &temp[0], &out[0], d-1);
	sub_fft(&in[h], &temp[h], &out[h], d-1);
	int p = 0;
	for (int i=0; i<h; i++)
	  { complex wkt = circle[p] * temp[h+i];
	    out[i] = temp[i] + wkt;
	    out[h+i] = temp[i] - wkt;
	    p += t;
	  }
      }
    else out[0] = in[0];
  }