Toby Smith

Goertzel Algorithm

9th September 2013

I’ve always wondered how complicated it would be to write the code that listens to the numbers you dial on your phone, and works out what numbers you pressed.
I often enter my credit card details over the phone, and think: “If someone had a good enough microphone to record the sounds, they could probably get my information!”

The fancy name for this is: Dual-tone multi-frequency signalling. When you press the buttons, two sine waves are generated with different frequencies. The receiver then has to measure these two frequencies and can work out what number you’ve pressed. I do love the similarity of this with formants (the way we create the different vowels sounds).

Whilst an FFT would easily do the trick, after a bit of research, it’s obvious that the most common algorithm to use is the Goertzel algorithm.
My time pocket principals say I should have just done an FFT and been done with it. The point here wasn’t to get a working decoder, but to understand how it is done, and to actually write some code that I can have.

I decided to work in scilab. I was tempted by C/C++, but I fancied learning a bit more about scilab (and not fiddling with audio libraries).

I created a couple of simple functions to generate some pure “tones” to test the algorithm, as well as adding noise, and gaps.

There’s plenty of examples of the Goertzel code out there. I started off doing different overlapping partitions to sample the audio, but in the end took a far simpler approach.
I’m not exactly satisfied with the speed, or with the way I’ve implemented the divisions. Ironically the first time I wrote it was probably a better method, but the second is definitely faster. So I certainly made some mistakes. Creating a matrix for all the frequencies at all the sample points probably was a bit silly… but I ran out of time.

In addition to the Goertzel I looked at dealing with noise, silence and series of numbers. The algorithm will sample the frequencies (by default) every 0.05 seconds, and will adjust to the sample rate of the audio data. If there isn’t a strong signal, then it will assume there’s no number being pressed. If the value is the same as the previous check, no new number is written out.

I tested it by recording my pressing buttons on a phone, and it actually worked out great! I was quite surprised.
All in all a good way to spend a few hours!

SciLab code here:
(I would imagine it would be fairly easy to port to Matlab).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89

// This creates an array of the goertzel's response.
// Use: data x, at frequency f, sampled at points defined by at.
function [res] = goertzelArr(x, f, sampleRate, at)
    a=0;
    b=0;
    c=0;
 
    val = 2 * cos(2 * %pi * f / sampleRate);
    res = zeros(1,length(at));
 
    nextAt = 1;
    nextTrig = at(1);
 
    for i = 1:length(x),
        a = val * b - c + x(i);
        c = b;
        b = a;
 
        if i == nextTrig then
            res(nextAt) = (b*b + c*c - b*c*val);            
            nextAt = nextAt + 1;
            b = 0;
            c = 0;
            if (nextAt > length(at) ) then
                return;
            end
            nextTrig = at(nextAt);
        end
    end
endfunction
 
// Main decoding function. 
// Call with: evalTone(AUDIO_VECTOR, SAMPLE_RATE_OF_AUDIO [, time between samples])
function [t] = evalToneAt(x,sampleRate,timebetween)
    // Default time between checks is 0.05
    [out, inp] = argn(0);
    if inp < 3 then, timebetween = 0.05, end
 
    // Work out how many samples per check, and build array    
    spaceSamples = timebetween * sampleRate;
    at = [spaceSamples:spaceSamples:length(x)];    
 
    // Relevant frequencies for number
    checkComb = [1209, 1336, 1477, 1633,697, 770, 852, 941];
    chars = ['1','2','3','A','4','5','6','B','7','8','9','C','*','0','#','D']; 
 
    // Create results array
    r = zeros(8,length(at));
 
    // Get results from goertzels
    for i=1:8
        r(i,:) = goertzelArr(x,checkComb(i),sampleRate,at);
    end
 
    lastVal = 0;
    t = '';
    // Loop through measurements
    for i=1:length(at)
        // Find best tone pair
        col = 1;
        row = 1;
        colval = r(1,i);
        rowval = r(5,i);
        for j=2:4
            if (r(j,i) > colval) then
                colval = r(j,i);
                col = j;
            end
            if (r(j+4,i) > rowval) then
                rowval = r(j+4,i);
                row = j;
            end
        end
        // This is the best candidate tone pair
        v = col + (row - 1) * 4;
 
        // If the magnitude doesn't satisfy contraints ignore.
        if (colval < spaceSamples | rowval < spaceSamples) then
            lastVal = 0;
            continue;
        end
 
        // If this is same tone as before. Ignore
        if (v ~= lastVal) then
            t = strcat([t,chars(v)]);
            lastVal = v;
        end
    end
endfunction

These functions create the sample tones:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

// Assistant function to create the tones
function [x] = soundChar(c, len, sampleRate)
 
    [out, inp] = argn(0);
    if inp < 3 then, sampleRate = 22050, end
    if inp < 2 then, len = 0.5, end    
 
    checkCol = [1209, 1336, 1477, 1633];
    checkRow = [697, 770, 852, 941];
 
    chars = ['1','2','3','A';'4','5','6','B';'7','8','9','C';'*','0','#','D'];
 
    sig1 = 0;
    sig2 = 0;
 
    for i = 1:size(chars,1)
        for j = 1:size(chars,2)
            if (chars(i,j) == c) then
                sig2 = checkRow(i);
                sig1 = checkCol(j);
            end
        end
    end
 
    v1 = len * 2 * %pi;
    v2 = len * sampleRate;
 
    s1 = sin(linspace(0, v1*sig1, v2));
    s2 = sin(linspace(0, v1*sig2, v2));
 
    x = s1 + s2;
 
endfunction
 
// Create a tone sequence eg: genTone('0123456789#*ABCD',0.2,44100,0.05);
// Will return an audio vector with each tone playing for 0.2 seconds, 
// with gaps in between of 0.05 seconds at a sample rate of 44100
function [x] = genTone(s, len, sampleRate, gap) 
    [out, inp] = argn(0);
    if inp < 4 then, gap = 0.2, end
    if inp < 3 then, sampleRate = 22050, end
    if inp < 2 then, len = 0.5, end   
 
    x = [];
    for i = 1:length(s)
        x = [x,soundChar(part(s,i), len, sampleRate),zeros(1,round(gap*sampleRate))];
    end
 
endfunction

Leave a Comment

http://www.thjsmith.com/feed">RSS Feed
  • Terms of use
  • Privacy Policy