```#!/usr/bin/python2

import numpy
import scipy.signal

tau = numpy.pi * 2
max_samples = 1000000
debug = False

# determine the clock frequency
# input: magnitude spectrum of clock signal (numpy array)
# output: FFT bin number of clock frequency
def find_clock_frequency(spectrum):
maxima = scipy.signal.argrelextrema(spectrum, numpy.greater_equal)
while maxima < 2:
maxima = maxima[1:]
if maxima.any():
threshold = max(spectrum[2:-1])*0.8
indices_above_threshold = numpy.argwhere(spectrum[maxima] > threshold)
return maxima[indices_above_threshold]
else:
return 0

def midpoint(a):
mean_a = numpy.mean(a)
high = numpy.ma.median(mean_a_greater)
low = numpy.ma.median(mean_a_less_or_equal)
return (high + low) / 2

# whole packet clock recovery
# input: real valued NRZ-like waveform (array, tuple, or list)
#        must have at least 2 samples per symbol
#        must have at least 2 symbol transitions
# output: list of symbols
def wpcr(a):
if len(a) < 4:
return []
b = (a > midpoint(a)) * 1.0
d = numpy.diff(b)**2
if len(numpy.argwhere(d > 0)) < 2:
return []
f = scipy.fft(d, len(a))
p = find_clock_frequency(abs(f))
if p == 0:
return []
cycles_per_sample = (p*1.0)/len(f)
clock_phase = 0.5 + numpy.angle(f[p])/(tau)
if clock_phase <= 0.5:
clock_phase += 1
symbols = []
for i in range(len(a)):
if clock_phase >= 1:
clock_phase -= 1
symbols.append(a[i])
clock_phase += cycles_per_sample
if debug:
print("peak frequency index: %d / %d" % (p, len(f)))
print("samples per symbol: %f" % (1.0/cycles_per_sample))
print("clock cycles per sample: %f" % (cycles_per_sample))
print("clock phase in cycles between 1st and 2nd samples: %f" % (clock_phase))
print("clock phase in cycles at 1st sample: %f" % (clock_phase - cycles_per_sample/2))
print("symbol count: %d" % (len(symbols)))
return symbols

# convert soft symbols into bits (assuming binary symbols)
def slice_bits(symbols):
symbols_average = numpy.average(symbols)
bits = (symbols >= symbols_average)
return numpy.array(bits, dtype=numpy.uint8)