Extending Goertzel algorithm to 24 kHz, 32 kHz and 48 kHz

I'm learning to implement Goertzel's algorithm to detect DTMF tones from recorded wave files. I got one implemented in python from here. It supports audio sampled at 8 kHz and 16 kHz. I would like to extend it to support audio files sampled at 24 kHz, 32 kHz and 48 kHz.

From the code I got from the link above, I see that the author has set the following precondition parameters/constants:

self.MAX_BINS = 8
if pfreq == 16000:
self.GOERTZEL_N = 210
self.SAMPLING_RATE = 16000
else:
self.GOERTZEL_N = 92
self.SAMPLING_RATE = 8000


According to this article, before one can do the actual Goertzel, two of the preliminary calculations are:

1. Decide on the sampling rate.
2. Choose the block size, N

So, the author has clearly set block size as 210 for 16k sampled inputs and 92 for 8k sampled inputs. Now, I would like to understand:

1. how the author has arrived at this block size?
2. what would be the block size for 24k, 32k and 48k samples?

Choice of $$N$$

Dual-tone multi-frequency (DTMF) tones should be of certain minimum duration (40 ms) and quality to be detected (reviewed in ITU TABLE A-1/Q.24). The detection algorithm that you cite uses a block duration of $$t=N/f_s,$$ where $$N$$ is the block length in samples and $$f_s$$ is the sampling frequency. The duration is either 210/(16 kHz) = 13.125 ms or 92/(8 kHz) = 11.5 ms. To ensure that the algorithm functions as expected, the block duration should be kept the same. I don't know why the lengths differ for the different sampling frequencies in the cited algo, but it should be okay to calculate $$N$$ as:

$$N = t\times f_s = 11.5\text{ ms}\times f_s,\tag{1}$$

which should give the same timing and frequency resolution as with the 8 kHz sampling frequency used in the original sources on Wikipedia by an unknown author.

Other changes needed

There are some worrisome large hard-coded constants (4.0e5 and 1.0e9) in the source code. Such hint to a lack of normalization:

#Check for minimum energy
if self.r[row] < 4.0e5:
msg = "energy not enough"
elif self.r[col] < 4.0e5:
msg = "energy not enough"
else:


#AT&T states that the noise must be 16dB down from the signal.
# Here we count the number of signals above the threshold and
#there ought to be only two.
if self.r[max_index] > 1.0e9:
t = self.r[max_index] * 0.158
else:
t = self.r[max_index] * 0.010


They may require changing if you change $$N,$$ by the same proportion, I think.

Fine-tuning $$N$$

Even if the Goertzel algorithm is not restricted to frequencies of a discrete Fourier transform (DFT) of length $$N,$$ (see this answer) there may be reason for tuning $$N$$ that is not explained by duration specifications. For sinusoidal input of frequency $$f$$, the magnitude of the same frequency detected by a Goertzel algorithm fluctuates slightly as function of time depending on the combination of frequency $$f$$ and $$N$$. For simplicity of the expressions that follow, frequency is expressed in radians as:

$$\omega=\frac{2\pi f}{f_s}\tag{2}$$

For initial phase $$\alpha$$ at the start of the window, the unnormalized detected magnitude will be:

$$A(\alpha) = \left|\sum_{k=0}^{N-1} e^{i \omega k}\cos(\omega k + \alpha)\right|\\ = \left|\frac{Ne^{-i\alpha}}{2} + \frac{e^{i(\alpha - \omega+\pi/2)} - e^{i(\alpha + 2N\omega - \omega + \pi/2)}}{4\sin(\omega)}\right|\tag{3}$$

The fluctuation can be seen in this example:

Figure 1. Blue: detected magnitude $$A$$ as function of initial phase $$\alpha$$ for sinusoidal input, with $$f = 1209\text{ Hz},$$ $$f_s = 8 kHz,$$ $$N=92$$. Red: extrema of $$A(\alpha)$$ have been marked.

An interesting observation is that there will no fluctuation of the Goertzel detected magnitude $$A$$ if the frequency $$\omega$$ is a multiple of $$\pi/N.$$ In contrast, DFT bins are at multiples of $$2\pi/N.$$ (Zaplata & Kasal 2005 discuss the fluctuation phenomenon, but do not recognize this difference between the Goertzel algorithm and DFT.)

For other frequencies, the $$\pi\text{-}$$periodic function $$A(\alpha)$$ gets its smallest and largest values (not necessarily in that order) at:

$$\alpha_1 = \frac{\omega(1 - N)}{2}\quad\text{and}\quad\alpha_2 = \frac{\omega(1 - N)}{2} + \frac{\pi}{2}\tag{4}$$

You could tune the integer $$N$$ to its nearby value which minimizes the maximum of $$\left|A(\alpha_1)-A(\alpha_2)\right|$$ over the set of DTMF frequencies. Actually, for $$f_s = 8\text{ kHz,}$$ a better choice than $$N=92$$ would be $$N=93:$$

Figure 2. Maximum peak-to-peak fluctuation $$\max\left|A(\alpha_1)-A(\alpha_2)\right|$$ over DTMF frequencies has a minimum at $$N=93$$ for $$f_s=8000\text{ Hz.}$$

A reason why $$N=92$$ might have been preferred is that in some older version they only allowed DFT bin frequencies. In that case $$N=93$$ would be a pretty bad choice, as some DTMF frequencies would fall almost just between successive bins, which is perfectly fine for the Goertzel algorithm (or generalized DFT with a frequency shift of $$1/2$$ binwidth), but almost the worst-case scenario when using DFT:

Figure 3. Maximum absolute frequency error over DTMF frequencies when detecting using DFT of length $$N$$, for $$f_s=8000\text{ Hz}$$ and expressed in units of DFT bin width.

The final word would be to test the performance in a test suite and to tune $$N$$ based on that.

Octave source code for Fig. 3

graphics_toolkit("gnuplot");
fs = 8000; #Sampling frequency
f = [1209, 1336, 1477, 1633, 697, 770, 852, 941]; #DTMF frequencies (same unit as fs)
N = 80:100; #DFT lengths analyzed
omega = 2*pi*f/fs;
y=omega'/(2*pi).*N;
z=round(y)-y;
plot(N, max(abs(z)), "x");
ylim([0, 0.5]);
xlabel("N");
ylabel('max|frequency error| (bin widths)');


References

Zaplata F, Kasal M, "Efficient Spectral Power Estimation on an Arbitrary Frequency Scale", Radioengineering, 2015, vol. 24: p. 178-184. DOI: 10.13164/re.2015.0178

• Great. That makes a lot of sense now. But how does one choose between block durations 13.125 ms and 11.5 ms? Commented Jan 2, 2019 at 9:33
• I think you meant N = t x fs or N = 13.125 ms x fs? Commented Jan 2, 2019 at 9:36
• Hey I'm John, the author. I wrote this back in college when I needed to know DTMF on what is now a very old Nokia phone. Back then the Wikipedia entry on the Goertzel algorithm had a C implementation that I copied, including all the constants. I don't have a strong math background, so I don't have any special knowledge on scaling up the sampling rate. The reasoning in Olli's answer makes sense to me, so I'd go with that. Thanks, and sorry I'm not more helpful. Commented Jan 2, 2019 at 16:03
• @etherton Your program and Olli's answer did help a lot. I've made it work with a minor modification in your code. Well, I only had a 48k sample and it works fine. I'm posting it in an answer post below. But I would still like to understand the concepts behind the block duration and other constants that Olli has pointed out. Commented Jan 3, 2019 at 3:55
• @OlliNiemitalo Thanks for the cues in your post. I have posted my observations on a trial and error basis. I would like to understand more on the underlying concept. I will dig into this and see if I can post a definite answer or wait till someone posts a more convincing answer. Commented Jan 3, 2019 at 4:24

With cues from Olli's answer, I have made the modifications in the program here.

The modification I did was that I kept the block duration as :

$$N = t\times f_s = 13.125\text{ ms}\times f_s,$$

So, according to this :

for $$f_s = 24\text{ kHz},$$ $$N = 315$$,

for $$f_s = 32\text{ kHz},$$ $$N = 420$$,

for $$f_s = 48\text{ kHz},$$ $$N = 630$$.

Which in the python code is translated as follows:

if pfreq == 48000:
# 48kHz samples
self.GOERTZEL_N = 630
self.SAMPLING_RATE = 48000
elif pfreq == 32000:
# 32kHz samples
self.GOERTZEL_N = 420
self.SAMPLING_RATE = 32000
elif pfreq == 24000:
# 24kHz samples
self.GOERTZEL_N = 315
self.SAMPLING_RATE = 24000
elif pfreq == 16000:
# 16kHz samples
self.GOERTZEL_N = 210
self.SAMPLING_RATE = 16000
else:
# 8kHz samples (default)
self.GOERTZEL_N = 92
self.SAMPLING_RATE = 8000


Though I have made it work, I'm not accepting this as an answer yet because I still want to know how the block duration is arrived at. Because, it seems to work with :

$$N = t\times f_s = 11.5\text{ ms}\times f_s,$$

as well. According to which :

for $$f_s = 16\text{ kHz},$$ $$N = 184$$,

for $$f_s = 24\text{ kHz},$$ $$N = 276$$,

for $$f_s = 32\text{ kHz},$$ $$N = 368$$,

for $$f_s = 48\text{ kHz},$$ $$N = 552$$.

I am waiting for an elaborate answer post which would help me understand the underlying concepts.

• Can you plot the detected tone level against the tone? I think you will find that the detection takes a little longer with longer block sizes, and that stays high a little longer after the original tone stops. And, the longer the block size, the "pickier" it will become - a longer block equates to a narrower filter. So, when deciding the block size, you have to consider sampling rate, minimum length of tones, and how close the tones will be (in terms of frequency.) You want it short to catch short bursts, but long enough to reliably separate your closest tones. Compromise.
– JRE
Commented Jan 3, 2019 at 12:09
• Also, rather than have a bunch of if/elses, I'd probably calculate the block length. Some constant like "13 milliseconds" then calculate the block length from that and the sampling rate.
– JRE
Commented Jan 3, 2019 at 12:11
• That's very interesting and enlightening. Thanks. I'll try plotting the tone and play around with N a bit. Commented Jan 4, 2019 at 9:35