first commit
This commit is contained in:
+211
@@ -0,0 +1,211 @@
|
||||
# Run with "ipython --matplotlib=qt receiver.py <file>.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])
|
||||
Reference in New Issue
Block a user