speed_from_tv.py

## Determine speed from frequency output of HB100 module
# Input file is logged with NI Datalogger

import sys
import numpy as np
from matplotlib import pyplot
from matplotlib.mlab import find

def speedFromTV(t, v):
	# t = time array / list
	# v = voltage array / list
	speedCoeff = 19.49 # see HB100 sensor application note

	dt = t[1] - t[0]
	windowHalfWidth_samples = 10 # 

	print "Analysis window length = %.2f milliseconds" % (2 * windowHalfWidth_samples * dt  * 1000)
	print "Theoretical highest speed %f" % (1 / (dt*2) / speedCoeff)

	speedInKmh_meanf = np.zeros(len(t))
	speedInKmh_medianf = np.zeros(len(t))

	for ti in xrange(windowHalfWidth_samples+1, len(v)-windowHalfWidth_samples):
		curBlock = np.array(v[ti - windowHalfWidth_samples:ti + windowHalfWidth_samples])
		waveLengths = 2*np.diff(find(np.abs(np.diff(curBlock)) > 1)) * dt
		if len(waveLengths) < 1:
			speedInKmh_medianf[ti] = 0.0
			speedInKmh_meanf[ti] = 0.0
		else:
			medianFreq = np.median(1.0 / waveLengths)
			meanFreq = np.mean(1.0 / waveLengths)
			speedInKmh_medianf[ti] = medianFreq / speedCoeff
			speedInKmh_meanf[ti] = meanFreq / speedCoeff


	return speedInKmh_meanf, speedInKmh_medianf


if __name__ == "__main__":
	raw = open(sys.argv[1]).read()

	t = np.array([float(row.split('\t')[0]) for row in raw.split('\r')[40:-1]])
	v = np.array([float(row.split('\t')[1]) for row in raw.split('\r')[40:-1]])

	speed_mean, speed_median = speedFromTV(t, v)

	pyplot.plot(t, speed_mean)
	pyplot.xlabel('t (s)')
	pyplot.ylabel('v (km/h)')
	pyplot.grid('on')

	pyplot.show()