Goertzel algorithm DTMF Detection

Goertzel algorithm is well known in the DSP domain. It is used to detect the presence of a frequency or a number of frequencies. This algorithm is popular as the calculation required for implementation is lesser as compared to other techniques such as DFT (FFT), when a smaller number of frequencies are to be detected.

Some projects that I have seen implementing Goertzel algorithm are: -

1. DTMF detection.

2. Tuning of guitar string

3. Detection of specific frequencies in audio signals.

My project uses PIC18F4520 and Goertzel algorithm to detect DTMF frequencies. The code is written entirely in "C". Even though, the implementation of my project is good, I do not guarantee that it can be used in a commercial application.

Understanding the Goertzel algorithm: -

The understanding of Goertzel algorithm can be made using the unit-circle diagram on z-plane (as shown below).

Usually for DTMF detection, the sample frequency is chosen to be 8000Hz. But I choose to use the sampling frequency to be 5908Hz. The reason for choosing this particular value will be explained later.

If we assume a system that has 2 conjugate poles (conjugate poles because, that effectively removes the complex term and make implementation possible using real world micro controller) on the unit circle, the output of such a system has a sharp peak at that frequency.

Sample frequency of 5908Hz falls on 2π radian on unit circle of z-plane.

Also, the frequencies beyond π is more than 1/2 of sampling frequency and hence falls beyond then Nyquist rate and hence we will not consider them.

Now, we have the following frequencies on DTMF: -

697Hz

770Hz

852Hz

941Hz

1209Hz

1336Hz

1477Hz

In order to place 697Hz on z-plane, we will normalize this with the sampling frequency: -

5908Hz ---> 2π

so, 697Hz = (2π/5908) * 697 = 0.74126 radians

Following the similar conversion, the values for each frequencies are: -

697Hz = 0.74126 radians

770Hz = 0.81889 radians

852Hz = 0.9061 radians

941Hz = 1.00713 radians

1209Hz = 1.28577 radians

1336Hz = 1.4208 radians

1477Hz = 1.57079 radians

If we place 2 conjugate poles for a system at frequency of 697Hz, and 2 zeros on origin, the z-plane diagram will look like this: -

And, the transfer functions of the system will be: -

(z – 0)(z – 0)

----------------------

(z – e^j0.74126)(z – e^-j0.74126)

In order to generalize our frequency, we will replace the term 0.74126 with "θ": -

(z – 0)(z – 0)

----------------------

(z – e^jθ)(z – e^-jθ)

z²

----------------------

z² – ze^-jθ - ze^jθ + 1

Using Eulers formula and replacing

e^jθ = Cosθ + jSinθ

and

e^-jθ = Cosθ - jSinθ

the equation becomes as below: -

z²

----------------------

z² – z(Cosθ – jSinθ) – z(Cosθ + jSinθ) + 1

z²

----------------------

z² – zCosθ – zjSinθ – zCosθ + zjSinθ + 1

z² Output

---------------------- = -------

z² – 2zCosθ + 1 Input

Multiplying numerator and denominator with "z^-2": -

z² * z^-2 Output

---------------------- = -------

(z² – 2zCosθ + 1) * z^-2 Input

1 Output

---------------------- = -------

1 – 2z^-1Cosθ + z^-2 Input

Replacing "z^-1" as one delay and "z^-2" as two delay: -

Input = Output - 2 * Cosθ * Prev_Output + Prev_Prev_Output

Output = Input + 2 * Cosθ * Prev_Output - Prev_Prev_Output

Above equation is the implementation of Goertzel algorithm. The term (2 * Cosθ) is called "coefficient" is pre-calculated for each frequencies: -

697Hz = 0.74126 radians = 1.4752

770Hz = 0.81889 radians = 1.3660

852Hz = 0.9061 radians = 1.2336

941Hz = 1.00713 radians = 1.0685

1209Hz = 1.28577 radians = 0.5623

1336Hz = 1.4208 radians = 0.2988

1477Hz = 1.57079 radians = 0

As you can see that the coefficient for "1477Hz" is zero and that basically save lot of calculation during detection of 1477Hz. This the reason why I choose a sampling frequency of 5908Hz.

These coefficients are placed as "#define" constants in "DTMF.h" file of my code. In order to speed up the calculation on PIC, I avoided any float calculation. Hence, the coefficients are multiplied by 256 to convert them into a integer. Later, I divide the result by 256 to get the actual value. This way, my calculations are way faster than float implementation.

The calculation for Goertzel is performed on a number of acquired ADC sample. I have chosen a sample number of 100 to perform the calculations. This is a balance between detection resolution of frequencies and processing time.

Once the calculations are performed on all the 100 samples, the magnitude of detection peak is calculated as: -

Magnitude = Prev_Prev_Output ² + Prev_Output² – (Coefficient * Prev_Prev_Output * Prev_Output)

If this magnitude is > than a threshold, then the frequency is present. Otherwise, the frequency is absent.

if(magnitude > THRESHOLD)

{

// Frequency is present

}

Understanding the Goertzel algorithm through excel sheet: -

For a more practical simulation of the algorithm, I have created an excel sheet. You can download the excel sheet and try different frequencies. Even though, the algorithm is used for DTMF detection, it can be used to detect any frequency provided the coefficients and sampling rates are correct.

Basic tests: -

I started my project by sampling a sine wave using PIC18F4520 ADC (channel - 0) triggered by Timer-3 (to make a precise time sampling rate). The Timer triggering ADC is absolutely essential as any jitter in sampling rate will degrade the performance of frequency detection. For testing purpose, I used a sample rate of 8KHz.

Screen shot for sample of 1KHz signal (@8KHz): -

Screen shot for sample of 400Hz signal (@8KHz): -

Screen shot for sample of 200Hz signal (@8KHz): -

In each of the above screen shot, count the number of points for 1 cycle should give the relation between (8000Hz/signal freq).

Design of hardware: -

Telephone line usually contains high voltages, specially during ringing, which can give some nasty shock. So isolating the telephone line from my circuit was essential. Hence, I had ordered the audio transformer (TY146P) from TRIAD. It is a audio transformer 1:1 with 600ohms impedance, which is perfect for telephone line coupling.

Circuit schematic: -

Explanation of hardware: -

1. 2 X 0.47uF X 450V capacitors are connected in parallel so that they block all DC voltage going into transformer. But they allow AC voltage like audio and DTMF signals to flow into the transformer.

2. 2K resistor is placed to reduce loading of telephone line.

3. TY-146P is a 600-600 ohms isolation transformer used for audio signal coupling.

4. 2 X 1N4007 diodes connected in antiparallel to limit the voltage that goes into the opamp.

5. LM358 -- one opamp used as voltage follower to generate 2.5V. This 2.5V is exactly the half of ADC reference used by PIC18F4520. Any AC signal to be measured by PIC18F4520 has to be DC offset by known voltage so that the negative sides of AC does not go below the Vss voltage (ground). In the software, 1/2 of ADC full range (in my case 129) is subtracted to remove the DC offset addition.

6. LM358 -- second opamp used as summer to add the 2.5V and input signal.

7. PIC18F4520 -- uses a 10MHz XTAL (internally X4PLL to generate 40MHz). Also, a LCD is interfaced to display the detected frequencies.

Once I rigged up the circuit, I did some measurement by connecting the circuit to my telephone line: -

FFT for signal while pressing "key - 1" (Freq = 697Hz, 1209Hz)

FFT for signal while pressing "key - 2" (Freq = 697Hz, 1336Hz)

FFT for signal while pressing "key - 3" (Freq = 697Hz, 1477Hz)

FFT for signal while pressing "key - 7" (Freq = 852Hz, 1209Hz)

Waveform when phone rings (without the protection diodes to limit volt level)

Please notice that when the phone rings, the voltage level goes very high. The above signal is captured after the isolation transformer, which had the capacitor and 2K resistor connected in primary side. The measured values is about +35V/-35V. If we do not use protection diode to clip the signal level, this high value of signal will certainly damage the analog input of microcontroller or the opamp.

Waveform when phone rings (with the antiparallel diodes 1N4007)

As we can see that the voltage level is safely limited to diode forward voltage of about +/-0.7V. The usage of diode does introduce some harmonics during the measurement of DTMF signal, but if the signal is kept low enough (within +/-0.7V), by increasing the value of input resistor R1 (from 2K to higher value), then the harmonics can be eliminated. Later, the signal can be amplified by amplifier to get sufficient level (if required). But if we use Goertzel algorithm, it will look for frequencies in specific range. Even if harmonics exists within the DTMF band, it will be ignored.

Advantages of using Goertzel/ software implementation of DTMF detection vs. dedicated chips: -

1. No extra hardware needed like encoder and decoder chips (like MT8870 and TP5089) and crystals. Using software implementation saves all that.

2. The equation for Goertzel can also be run from a timer (reverse Goertzel) to generate sinusoidal frequency that acts like DTMF generator. Hence, the encoder can also be implemented in software.

3. If a hardware design needs isolation, then the isolation transformer is required both with hardware implementation or Goertzel implementation. Same is for high voltage capacitors.

4. The concept for Goertzel can be used to detect any frequency which is not possible with hardware chips, as they are built to detect DTMF frequencies only.

5. Software implementation gives a fair idea about DSP concepts -- this is the main reason for me to like the software implementation.

Limitations of my project: -

1. Because of processing power limits, I am not detecting the second harmonics for DTMF. This is usually done in commercial implementation so as to eliminate any chance of false detection due to music or speech. Music or speech cannot generate pure sinusoids and hence contains second harmonics. But DTMF signals are pure sinusoids and hence does not contain second harmonics.

2. The peak detection might detect a false peak as there exist some ripples near the real peak. These ripples can be eliminated in a commercial implementation by using window function. I did not use window function due to processing power limits.

MIPS requirement: -

PIC18F4520 runs at a max frequency of 40MHz and internally, the clock is divided by 4 to get 10MIPS.

The implementation of Goertzel acquires a buffer of 100sample, sampled at a rate of 5908Hz.

Time to acquire 100 samples = 100 * (1/5908) = 16.92millisec

Time to perform Goertzel decoding (measured by toggling a GPIO on oscilloscope) = 9.40millisec

MIPS required = (9.40/16.92) * 10 = 5.55MIPS.

This means that my implementation of Goertzel for DTMF decoding takes > 50% of processing. More optimization can be done in the code (if I spend more time, which I cannot) or by converting C code into assembly code (which I do not like unless forced upon).

Final application: -

I wanted to create a home automation system and initially thought of using xigbee modules. But they are expensive. Later, it came to my mind that my house has a telephone network which goes to all the rooms. So I decided to build low cost DTMF boards to do home automation. The information for home automation will be sent as DTMF codes and this is the end use of my project.

Demonstration: -