How to detect complex signals, e.g. DTMF tones

Part One: Using the Cross Correlation

In many cases you come across a mixture of signals, and you want to check if the sequence of signals contains one or more elements of interest. Obviously, you cannot use a frequency counter to detect a certain signal in a noisy environment. And sometimes a Fourier analysis is not handy and even heavily time consuming. So a cross-correlation with patterns of what you are looking for can be helpful.

In the example described below we apply DTMF tones produced by a tone remote controller for answering machines (Fig. 1) to a microphone (we took the one equipped with the MAX9814 amplifier, Fig. 2) connected to an Arduino in order do find out which button was pressed.

For the younger readers: devices like this were used once upon a time to check your tape-driven answering machine from an elderly telephone booth still equipped with the old-fashioned rotary phone dial. Nowadays, you do not need a separate device like this any longer as nearly all phones offer to produce the DTMF tones.

To find the level of the sine tone you are looking for this line will give you the answer:

for (int i = 0; i < N; i++) 
   c = c + x[i] * sin(TWO_PI * dt / N * f[j] * i);

where x[i] are the sampled data and f[j] is the list of the frequencies you are looking for.

(In math lessons, you will be taught that the sine functions build an orthogonal system of functions.)

(To decode DTMF with software, often the Goertzel algorithm is used which is very fast and effective, but less easy to implement and explain.)

When you invoke the Serial plotter, a window like this will show up:

In this case, button “1” was pressed which produces 697 Hz and 1209 Hz.

Actually, most of the code was written to enable the Serial plotter to produce this colourful kind of diagram.

/*
  connect microphone with Max9814 to pins A0-A4
  use remote control for answering machine
  or mobile phone and press some buttons
  the SERIAL PLOTTER at its best!
*/

const byte input = A1;
const byte Vcc = A3;
const byte GND = A4;

void setup() {
  Serial.begin(115200);
  Serial.println(__FILE__);
  pinMode(Vcc, OUTPUT);
  pinMode(GND, OUTPUT);
  digitalWrite(Vcc, HIGH);
}

void loop() {
  const int N = 500;
  int x[N];
  long t1 = micros();
  for (int i = 0; i < N; i++) x[i] = analogRead(input);
  long t2 = micros();
  float dt = (t2 - t1) * 1E-6; // seconds
  long sum = 0;
  for (int i = 0; i < N; i++) sum = sum + x[i];
  int avg = sum / N;
  for (int i = 0; i < N; i++) x[i] = x[i] - avg; // normalize
  int f[] = { 697, 770, 852, 941, 1209, 1336, 1477, 1633 };
  byte n = (sizeof f) / 2;
  // gain some space for new values:
  Serial.println("0 0 0 0 0 0 0 10000"); // outsmart autoscaling
  Serial.println("0 0 0 0 0 0 0 0");
  Serial.println("0 0 0 0 0 0 0 0");
  float b = TWO_PI * dt / N; // for speed only
  for (byte j = 0; j < n; j++) {
    float d = b * f[j];
    float c = 0;
    for (int i = 0; i < N; i++)
      //================================+
      c = c + x[i] * sin(d * i);     // |
    // =================================+
    c = abs(c);
    for (byte i = 0; i < 58; i++) {
      for (byte k = 0; k < n; k++)
        if (k != j) Serial.print(" 0");
        else { // use separate colors
          Serial.print(" ");
          Serial.print(c);
        }
      Serial.println();
    }
    Serial.println("0 0 0 0 0 0 0 0");
    Serial.println("0 0 0 0 0 0 0 0");
    Serial.println("0 0 0 0 0 0 0 0");
    Serial.println("0 0 0 0 0 0 0 0");
  }
  delay(500);
}

Part Two: Fourier analysis - direct calculation compared to Fast Fourier Transform

Now let's compare the original algorithm with the one developed by Cooley and Tukey in 1965. The main program is the same for both of them. Use a jumper on pins 8 & 9 to select the FFT.

/* Fourier-Analysis
  For ARDUINO UNO
  or ARDUINO MEGA2560.
  Insert the Microphone module MAX9814
  into the pins A0-A4 (GND=A4).
  Optional: Jumper to d8-d9.
  Use the Serial Plotter for display.
*/
const byte audioIn = A1;
const byte Vcc = A3;
const byte GND = A4;
const byte selPin1 = 8;
const byte selPin2 = 9;
#if defined(__AVR_ATmega2560__)
const byte BITS = 8;
#else
// UNO does not offer enough RAM memory:
const byte BITS = 7;
#endif
const int N = 1 << BITS;
int s[N]; // samples of audioIn
float a[N], b[N]; // coefficients
#include "std.h"
#include "fft.h"

void setup() {
  Serial.begin(115200);
  // When using ARDUINO UNO all values
  // will be printed twice to use
  // the plotter window.
  pinMode(GND, OUTPUT);
  pinMode(Vcc, OUTPUT);
  digitalWrite(Vcc, HIGH);
  pinMode(selPin1, OUTPUT); // set to GND
  pinMode(selPin2, INPUT_PULLUP);
  // Wait until power to the module 
  // is stable:
  delay(65);
}

void loop() {
  // read samples :
  for (int k = 0; k < N; k++)
    s[k] = analogRead(audioIn);
  long mw = 0; // mean
  for (int k = 0; k < N; k++)
    mw = mw + s[k];
  mw = mw / N;
  // print the samples (blue):
  for (int k = 0; k < N; k++) {
    Serial.print(s[k] - mw);
    Serial.println(" 0");
#if defined(__AVR_ATmega328P__)
    Serial.print(s[k] - mw);
    Serial.println(" 0");
#endif
  }
  // calculate the coefficients:
  /* --------------------------------- */
  // select the algorithm:
  if (digitalRead(selPin2) == HIGH)
    fourierStd(); // no jumper
  else
    fourierFft(); // jumper set
  /* --------------------------------- */
  float f = TWO_PI / N; // normalize
  // print abs (red):
  for (int m = 1; m < N; m++) {
    float v = sqrt(sq(a[m]) + sq(b[m])) * f;
    Serial.print("0 ");
    Serial.println(v);
#if defined(__AVR_ATmega328P__)
    Serial.print("0 ");
    Serial.println(v);
#endif
  }
  delay(1000);
}

Now this is the original version:

// filename should be: "std.h"
void fourierStd() {
  /*
    see Wikipedia:
    Jean-Baptiste Joseph Fourier was a
    French mathematician and physicist best
    known for initiating the investigation
    of Fourier series, which eventually
    developed into Fourier analysis
    and harmonic analysis.
    Here the direct calculation:
  */
  // calculate the coefficients:
  for (int m = 1; m < N; m++) {
    a[m] = 0;
    b[m] = 0;
    for (int k = 0; k < N; k++) {
      float arg = TWO_PI * m * k / N;
      a[m] = a[m] + s[k] * cos(arg);
      b[m] = b[m] + s[k] * sin(arg);
    }
  }
}

Quite short - isn't it? But it is O(n*n).

And this is the improved version:

// filename should be: "fft.h"
void fourierFft() {
  // Fast Fourier Transformation:
  // calculate FFT:
  // Algorithm by James William Cooley
  // uad John Tukey (IBM), 1965
  for (int k = 0; k < N; k++) {
    // copy int to float:
    a[k] = s[k];
    b[k] = 0;
  }
  float tmp;
  int j = 0;
  for (int i = 0; i < N - 1; i++) {
    if (i < j) {
      tmp = a[i];
      a[i] = a[j];
      a[j] = tmp;
      tmp = b[i];
      b[i] = b[j];
      b[j] = tmp;
    }
    // bit-reverse:
    int k = (N >> 1);
    while (k <= j) {
      j -= k;
      k >>= 1;
    }
    j += k;
  }
  float c1 = -1;
  float c2 =  0;
  int l2 = 1;
  for (int l = 0; l < BITS; l++) {
    int l1 = l2;
    l2 <<= 1;
    float u1 = 1;
    float u2 = 0;
    for (int j = 0; j < l1; j++) {
      for (int i = j; i < N; i += l2) {
        int i1 = i + l1;
        float t1 = u1 * a[i1] - u2 * b[i1];
        float t2 = u1 * b[i1] + u2 * a[i1];
        a[i1] = a[i] - t1;
        b[i1] = b[i] - t2;
        a[i] += t1;
        b[i] += t2;
      }
      tmp = u1 * c1 - u2 * c2;
      u2 = u1 * c2 + u2 * c1;
      u1 = tmp;
    }
    c2 = -sqrt((1 - c1) / 2);
    c1 =  sqrt((1 + c1) / 2);
  }
}

Much longer, but much faster: O(n * log(n) )

You can insert and remove the jumper whenever you want. Be aware that the original code will take very long when using the ATmega2560 as it supports twice the sample size.

And this is what you get:

The left half of the diagram shows the time series of the microphone signal (blue). The right half shows the result of the Fourier transform (red). Be aware that the Fourier algorithm lists the spectrum upwards from zero to the Nyquist frequency and then back to zero, so you get all the peaks twice. For calibration of the diagram you might use a tunig fork.

The Jumper

The program offers the selection of two different algorithms. This could be done by uploading different sketches which takes long time. Or by means of a switch. I decided to use a jumper which is very easy to handle. My jumper consists of a resistor of some 100 ohms. Why this? In class, students might attach the jumper to +Vcc and GND, and this could easily kill the USB port of the computer.

When using this kind of jumper nothing serious can happen.