diff --git a/TestSelcall.wav b/TestSelcall.wav new file mode 100644 index 0000000..8d7bbb5 Binary files /dev/null and b/TestSelcall.wav differ diff --git a/receiver.py b/receiver.py new file mode 100644 index 0000000..ee01dce --- /dev/null +++ b/receiver.py @@ -0,0 +1,211 @@ +# Run with "ipython --matplotlib=qt receiver.py .wav" +# +from __future__ import print_function +import sys +import numpy as np +from scipy import signal +from scipy.io.wavfile import read +from scipy.signal import butter, lfilter +from math import log10, floor + +from matplotlib import pyplot as plt +from mpl_toolkits.mplot3d import Axes3D + +FRAME_TIME = 0.04 # Frame time in seconds + +TONES = [1124, + 1197, + 1275, + 1358, + 1446, + 1540, + 1640, + 1747, + 1860, + 1981, + 2110] + +ALPHABET = ['1', + '2', + '3', + '4', + '5', + '6', + '7', + '8', + '9', + '0', + 'E'] + + +reffreq = {} +reffreq['1'] = 1124 +reffreq['2'] = 1197 +reffreq['3'] = 1275 +reffreq['4'] = 1358 +reffreq['5'] = 1446 +reffreq['6'] = 1540 +reffreq['7'] = 1640 +reffreq['8'] = 1747 +reffreq['9'] = 1860 +reffreq['0'] = 1981 +reffreq['e'] = 2110 + +FILTER_LEN = 1000 # Samples + + +# Shamelessly lifted from +# https://scipy.github.io/old-wiki/pages/Cookbook/ButterworthBandpass +def butter_bandpass(lowcut, highcut, fs, order=5): + nyq = 0.5 * fs + low = lowcut / nyq + high = highcut / nyq + b, a = butter(order, [low, high], btype='band') + return b, a + + +def butter_bandpass_filter(data, lowcut, highcut, fs, order=5): + b, a = butter_bandpass(lowcut, highcut, fs, order=order) + y = lfilter(b, a, data) + return y + + +# Tone synthesis +def note(freq, cycles, amp=32767.0, rate=44100): + len = cycles * (1.0/rate) + t = np.linspace(0, len, int(len * rate)) + if freq == 0: + data = np.zeros(int(len * rate)) + else: + data = np.sin(2 * np.pi * freq * t) * amp + return data.astype(int) + +def checktrain(train): + # Look for 5 definitive tones in a row + for idx,[tone,corr] in enumerate(train): + #print(idx,tone,corr) + ener = 0 + val = [] + + if(idx>4): + for i in range(idx-5,idx): + tone,corr = train[i] + val.append(tone) + ener += corr + + if(ener > 54): + print(ener," value:",val) + + + + return 0 + + +# analyze wav file by chunks +def receiver(file_name): + try: + sig_rate, sig_noise = read(file_name) + except Exception: + print('Error opening {}'.format(file_name)) + return + + print('file: ', file_name, ' rate: ', sig_rate, ' len: ', len(sig_noise)) + + if sig_rate == 44100: + decimate = 4 # rate = 11025, Fmax = 5512.5 Hz + elif sig_rate == 48000: + decimate = 5 # rate = 9600, Fmax = 4800 Hz + elif sig_rate == 22050: + decimate = 2 # rate = 11025, Fmax = 5512.5 Hz + elif sig_rate == 11025: + decimate = 1 # rate = 11025, Fmax = 5512.5 Hz + else: + print('Sample rate {} not supported.'.format(sig_rate)) + return + + if decimate > 1: + sig_noise = signal.decimate(sig_noise, decimate) + sig_rate = sig_rate / decimate + print('length after decimation: ', len(sig_noise)) + + frame_len = int(sig_rate * FRAME_TIME) + frames = int(floor((len(sig_noise) / frame_len) + 1)) + + sig_noise = butter_bandpass_filter(sig_noise, + 1000, + 2200, + sig_rate, + order=8) + + template = [] + for tone in range(0, len(TONES)): + template.append(note(TONES[tone], frame_len, rate=sig_rate)) + + # See http://stackoverflow.com/questions/23507217/ + # python-plotting-2d-data-on-to-3d-axes + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + + y = np.arange(len(TONES)) + print(' Index 1 2 3 4 5 6 7', end='') + print(' 8 9 0 E Avg') + + x = range(0, frames) + X, Y = np.meshgrid(y, x) + Z = np.zeros((len(x), len(y))) + + + hist = [] + for frame in range(0, frames): + + beg = frame * frame_len + end = (frame+1) * frame_len + + corr = np.zeros(len(TONES)) + + for tone in range(0, len(TONES)): + corr[tone] = log10(np.abs(signal.correlate(sig_noise[beg:end],template[tone],mode='same')).sum()) + Z[frame, tone] = corr[tone] + + max1 = 0.0 + # Find most likely tone in set + for tone in range(0, len(TONES)): + if corr[tone] > max1: + max1 = corr[tone] + max1idx = tone + + hist.append([max1idx,max1]) + + print('{0:6d}: '.format(frame), end='') + avg = np.mean(corr) + for tone in range(0, len(TONES)): + if tone == max1idx: + print('[{0:2.2f}]'.format(corr[tone]), end='') + else: + if corr[tone] > 1: + print(' {0:2.2f} '.format(corr[tone]), end='') + else: + print(' . ', end='') + + print(' {0:2.2f}'.format(avg)) + + checktrain(hist) + ax.plot_surface(X, Y, Z, rstride=1, cstride=1000, color='w', shade=True, + lw=.5) + # ax.plot_wireframe(X, Y, Z, rstride=1, cstride=1000, lw=.5) + + ax.set_title(file_name) + ax.set_xlabel("Tone") + ax.set_ylabel("Frame") + ax.set_zlabel("Log Correlation") + + ax.set_zlim(10.0, 15.0) + ax.set_ylim(0, frames) + + ax.view_init(30, -130) + + plt.show() + + +if __name__ == "__main__": + receiver(sys.argv[1])