From e7e99516612d912d7c6a04537eed3124eeb32070 Mon Sep 17 00:00:00 2001 From: "Frank J. T. Wojcik" Date: Thu, 16 May 2024 14:46:05 -0700 Subject: [PATCH 1/8] Whitespace fixes --- dtmf.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dtmf.py b/dtmf.py index 7c17e83..6437779 100755 --- a/dtmf.py +++ b/dtmf.py @@ -51,7 +51,7 @@ print ("Warning: The sound is not mono and not stereo ("+str(data.shape[1])+" canals)... so the -r option was ignored.") else: - if len(data.shape) == 2: + if len(data.shape) == 2: data = data.sum(axis=1) # stereo precision = args.i @@ -88,15 +88,15 @@ plt.plot(signal) plt.xticks([]) plt.yticks([]) - - + + frequencies = np.fft.fftfreq(signal.size, d=1/fps) amplitudes = np.fft.fft(signal) # Low i_min = np.where(frequencies > 0)[0][0] i_max = np.where(frequencies > 1050)[0][0] - + freq = frequencies[i_min:i_max] amp = abs(amplitudes.real[i_min:i_max]) From 67e17b241d5a89a2cf2e8f4045935d14780c488e Mon Sep 17 00:00:00 2001 From: "Frank J. T. Wojcik" Date: Thu, 16 May 2024 14:49:44 -0700 Subject: [PATCH 2/8] Append mirrored copy of signal for FFT analysis --- dtmf.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/dtmf.py b/dtmf.py index 6437779..71b9ce5 100755 --- a/dtmf.py +++ b/dtmf.py @@ -73,6 +73,10 @@ try: for i in range(0, len(data)-step, step): signal = data[i:i+step] + # Since Fourier transforms work best on periodic signals, we can + # fake it for this slice of signal by appending a mirrored + # copy. This helps limit the noise in low-frequency bins. + signal = np.append(signal, np.flip(signal)) if debug: plt.subplot(311) From d76c156c1bc58432f6f430719e9e034584fcd470 Mon Sep 17 00:00:00 2001 From: "Frank J. T. Wojcik" Date: Thu, 16 May 2024 14:55:38 -0700 Subject: [PATCH 3/8] Add option to slide window by less than full width A default of 3 is a good improvement in my tests --- dtmf.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/dtmf.py b/dtmf.py index 71b9ce5..0b8046f 100755 --- a/dtmf.py +++ b/dtmf.py @@ -15,6 +15,7 @@ parser.add_argument("-d", "--debug", help="show graphs to debug", action="store_true") parser.add_argument("-t", type=int, metavar="F", help="acceptable frequency error (in hertz, 20 by default)", default=20) parser.add_argument("-i", type=float, metavar='T', help="process by T seconds intervals (0.04 by default)", default=0.04) +parser.add_argument("-f", type=float, metavar="FRAC", help="process by shifting interval by 1/FRAC (3.0 by default)", default=3.0) parser.add_argument('file', type=argparse.FileType('r')) @@ -58,7 +59,9 @@ duration = len(data)/fps -step = int(len(data)//(duration//precision)) +window = int(len(data)//(duration//precision)) + +step = int(window//args.f) debug = args.debug verbose = args.verbose @@ -71,8 +74,8 @@ print ("0:00 ", end='', flush=True) try: - for i in range(0, len(data)-step, step): - signal = data[i:i+step] + for i in range(0, len(data)-window, step): + signal = data[i:i+window] # Since Fourier transforms work best on periodic signals, we can # fake it for this slice of signal by appending a mirrored # copy. This helps limit the noise in low-frequency bins. @@ -86,7 +89,7 @@ plt.xticks([]) plt.yticks([]) plt.axvline(x=i, linewidth=1, color='red') - plt.axvline(x=i+step, linewidth=1, color='red') + plt.axvline(x=i+window, linewidth=1, color='red') plt.subplot(312) plt.title("analysed frame") plt.plot(signal) @@ -154,9 +157,9 @@ plt.xlabel(txt) - t = int(i//step * precision) + t = int(i//window * precision) - if verbose and t > int((i-1)//step * precision): + if verbose and t > int((i-step)//window * precision): m = str(int(t//60)) s = str(t%60) s = "0"*(2-len(s)) + s From cf5b0e379861908868adcf9f51887fd5f719dfe1 Mon Sep 17 00:00:00 2001 From: "Frank J. T. Wojcik" Date: Thu, 16 May 2024 16:40:20 -0700 Subject: [PATCH 4/8] Reorder plt commands, to allow for subroutine creation --- dtmf.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/dtmf.py b/dtmf.py index 0b8046f..d9a25cf 100755 --- a/dtmf.py +++ b/dtmf.py @@ -95,6 +95,8 @@ plt.plot(signal) plt.xticks([]) plt.yticks([]) + plt.subplot(313) + plt.title("Fourier transform") frequencies = np.fft.fftfreq(signal.size, d=1/fps) @@ -118,8 +120,6 @@ best = f if debug: - plt.subplot(313) - plt.title("Fourier transform") plt.plot(freq, amp) plt.yticks([]) plt.annotate(str(int(lf))+"Hz", xy=(lf, max(amp))) @@ -150,6 +150,7 @@ hf = best if debug: + plt.yticks([]) if lf == 0 or hf == 0: txt = "Unknown dial tone" else: From 59c47a1fb1a853ffe0b65e23bf13c827f5f5cebf Mon Sep 17 00:00:00 2001 From: "Frank J. T. Wojcik" Date: Thu, 16 May 2024 16:47:08 -0700 Subject: [PATCH 5/8] Create findbest() subroutine --- dtmf.py | 50 +++++++++++++++++++------------------------------- 1 file changed, 19 insertions(+), 31 deletions(-) diff --git a/dtmf.py b/dtmf.py index d9a25cf..0f08bcb 100755 --- a/dtmf.py +++ b/dtmf.py @@ -73,6 +73,23 @@ if verbose: print ("0:00 ", end='', flush=True) +def findbest(freq, amp, maxdelta, freqlist): + maxf = freq[np.where(amp == max(amp))[0][0]] + + delta = maxdelta + best = 0 + + for f in freqlist: + if abs(maxf-f) < delta: + delta = abs(maxf-f) + best = f + + if debug: + plt.plot(freq, amp) + plt.annotate(str(int(maxf))+"Hz", xy=(maxf, max(amp))) + + return best + try: for i in range(0, len(data)-window, step): signal = data[i:i+window] @@ -109,22 +126,7 @@ freq = frequencies[i_min:i_max] amp = abs(amplitudes.real[i_min:i_max]) - lf = freq[np.where(amp == max(amp))[0][0]] - - delta = args.t - best = 0 - - for f in [697, 770, 852, 941]: - if abs(lf-f) < delta: - delta = abs(lf-f) - best = f - - if debug: - plt.plot(freq, amp) - plt.yticks([]) - plt.annotate(str(int(lf))+"Hz", xy=(lf, max(amp))) - - lf = best + lf = findbest(freq, amp, args.t, [697, 770, 852, 941]) # High i_min = np.where(frequencies > 1100)[0][0] @@ -133,21 +135,7 @@ freq = frequencies[i_min:i_max] amp = abs(amplitudes.real[i_min:i_max]) - hf = freq[np.where(amp == max(amp))[0][0]] - - delta = args.t - best = 0 - - for f in [1209, 1336, 1477, 1633]: - if abs(hf-f) < delta: - delta = abs(hf-f) - best = f - - if debug: - plt.plot(freq, amp) - plt.annotate(str(int(hf))+"Hz", xy=(hf, max(amp))) - - hf = best + hf = findbest(freq, amp, args.t, [1209, 1336, 1477, 1633]) if debug: plt.yticks([]) From 8e2150e3fbb875e70f5db78319504458dcaec340 Mon Sep 17 00:00:00 2001 From: "Frank J. T. Wojcik" Date: Thu, 16 May 2024 18:09:45 -0700 Subject: [PATCH 6/8] Improve freq matching algorithm --- dtmf.py | 71 ++++++++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 55 insertions(+), 16 deletions(-) diff --git a/dtmf.py b/dtmf.py index 0f08bcb..c977df9 100755 --- a/dtmf.py +++ b/dtmf.py @@ -16,6 +16,8 @@ parser.add_argument("-t", type=int, metavar="F", help="acceptable frequency error (in hertz, 20 by default)", default=20) parser.add_argument("-i", type=float, metavar='T', help="process by T seconds intervals (0.04 by default)", default=0.04) parser.add_argument("-f", type=float, metavar="FRAC", help="process by shifting interval by 1/FRAC (3.0 by default)", default=3.0) +parser.add_argument("-s", type=float, metavar="SD", help="ignore signals less than SD above median (3.0 by default)", default=3.0) +parser.add_argument("-m", type=float, metavar="MIN", help="ignore hf signals less than 1/MIN above lf (3.0 by default)", default=3.0) parser.add_argument('file', type=argparse.FileType('r')) @@ -73,22 +75,56 @@ if verbose: print ("0:00 ", end='', flush=True) -def findbest(freq, amp, maxdelta, freqlist): - maxf = freq[np.where(amp == max(amp))[0][0]] - - delta = maxdelta - best = 0 - - for f in freqlist: - if abs(maxf-f) < delta: - delta = abs(maxf-f) - best = f +def findbest(freq, amp, maxdelta, minamp, freqlist): + # Look for signals 3+ deviations (by default) above the median + medf = np.median(amp) + sdf = np.std(amp) + sdbound = medf + args.s * sdf + # Also obey any given minimum amplitude + bound = max(minamp, sdbound) + + # If there is no such signal, just try the strongest + fpos = np.where(amp >= bound)[0] + if len(fpos) == 0: + fpos = np.where(amp == max(amp))[0] + + # Find the strongest signal that is close enough (within maxdelta) of + # one we're looking for + bestamp = 0 + bestidealf = 0 + bestrealf = 0 + + for pos in fpos: + if bestamp >= amp[pos]: + continue + delta = maxdelta + realf = freq[pos] + for idealf in freqlist: + if abs(realf-idealf) < delta: + delta = abs(realf-idealf) + bestidealf = idealf + if delta < maxdelta: + bestamp = amp[pos] + bestrealf = realf + + # If we found no matching ideal signal, return the strongest signal given + if bestamp == 0: + bestamp = max(amp) if debug: - plt.plot(freq, amp) - plt.annotate(str(int(maxf))+"Hz", xy=(maxf, max(amp))) - - return best + plt.plot(freq, amp, color='blue' if minamp == 0 else 'green') + plt.plot([freq[0], freq[-1]], [sdbound, sdbound], + linestyle='dashed', color='orange') + if minamp: + plt.plot([freq[0], freq[-1]], [minamp, minamp], + linestyle='dashed', color='red') + for f in freq[fpos]: + a = amp[np.where(freq == f)] + plt.annotate(str(int(f))+"Hz", + xy=(f, a), + color='green' if f == bestrealf else 'orange') + + return [bestidealf, bestamp] try: for i in range(0, len(data)-window, step): @@ -126,7 +162,7 @@ def findbest(freq, amp, maxdelta, freqlist): freq = frequencies[i_min:i_max] amp = abs(amplitudes.real[i_min:i_max]) - lf = findbest(freq, amp, args.t, [697, 770, 852, 941]) + [lf, lfamp] = findbest(freq, amp, args.t, 0, [697, 770, 852, 941]) # High i_min = np.where(frequencies > 1100)[0][0] @@ -135,7 +171,10 @@ def findbest(freq, amp, maxdelta, freqlist): freq = frequencies[i_min:i_max] amp = abs(amplitudes.real[i_min:i_max]) - hf = findbest(freq, amp, args.t, [1209, 1336, 1477, 1633]) + minamp = lfamp // args.m; + + [hf, hfamp] = findbest(freq, amp, args.t, minamp, [1209, 1336, 1477, 1633]) + if debug: plt.yticks([]) From 10f047f34b2539b2c0312f01dd2eac0a1622ae7f Mon Sep 17 00:00:00 2001 From: "Frank J. T. Wojcik" Date: Fri, 17 May 2024 12:54:35 -0700 Subject: [PATCH 7/8] Gather frequencies from dtmf dict --- dtmf.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/dtmf.py b/dtmf.py index c977df9..1d8610e 100755 --- a/dtmf.py +++ b/dtmf.py @@ -69,6 +69,12 @@ verbose = args.verbose c = "" +lfset = set() +hfset = set() +for key in dtmf: + lfset.add(key[0]) + hfset.add(key[1]) + if debug: print("Warning:\nThe debug mode is very uncomfortable: you need to close each window to continue.\nFeel free to kill the process doing CTRL+C and then close the window.\n") @@ -162,7 +168,7 @@ def findbest(freq, amp, maxdelta, minamp, freqlist): freq = frequencies[i_min:i_max] amp = abs(amplitudes.real[i_min:i_max]) - [lf, lfamp] = findbest(freq, amp, args.t, 0, [697, 770, 852, 941]) + [lf, lfamp] = findbest(freq, amp, args.t, 0, list(lfset)) # High i_min = np.where(frequencies > 1100)[0][0] @@ -173,7 +179,7 @@ def findbest(freq, amp, maxdelta, minamp, freqlist): minamp = lfamp // args.m; - [hf, hfamp] = findbest(freq, amp, args.t, minamp, [1209, 1336, 1477, 1633]) + [hf, hfamp] = findbest(freq, amp, args.t, minamp, list(hfset)) if debug: From 50883c8c63eb5e8caea2d44a0055ac761bc8613c Mon Sep 17 00:00:00 2001 From: "Frank J. T. Wojcik" Date: Fri, 17 May 2024 13:07:37 -0700 Subject: [PATCH 8/8] Tweak default params for better performance --- dtmf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dtmf.py b/dtmf.py index 1d8610e..bef8dfb 100755 --- a/dtmf.py +++ b/dtmf.py @@ -14,8 +14,8 @@ parser.add_argument("-r", "--right", help="right channel only (if the sound is stereo)", action="store_true") parser.add_argument("-d", "--debug", help="show graphs to debug", action="store_true") parser.add_argument("-t", type=int, metavar="F", help="acceptable frequency error (in hertz, 20 by default)", default=20) -parser.add_argument("-i", type=float, metavar='T', help="process by T seconds intervals (0.04 by default)", default=0.04) -parser.add_argument("-f", type=float, metavar="FRAC", help="process by shifting interval by 1/FRAC (3.0 by default)", default=3.0) +parser.add_argument("-i", type=float, metavar='T', help="process by T seconds intervals (0.05 by default)", default=0.05) +parser.add_argument("-f", type=float, metavar="FRAC", help="process by shifting interval by 1/FRAC (2.0 by default)", default=2.0) parser.add_argument("-s", type=float, metavar="SD", help="ignore signals less than SD above median (3.0 by default)", default=3.0) parser.add_argument("-m", type=float, metavar="MIN", help="ignore hf signals less than 1/MIN above lf (3.0 by default)", default=3.0)