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
