From 4b267b673083a244e54dba8c90d5dcd0e510b61b Mon Sep 17 00:00:00 2001 From: Sebastien Bourdeauducq Date: Sun, 5 Jan 2020 19:31:07 +0800 Subject: [PATCH] cleanup, fix SDR buffer overflows --- stabilize.py | 115 ++++++++++++++++++++++---------------------------- stabilizer.py | 39 +++++++++++++++++ 2 files changed, 89 insertions(+), 65 deletions(-) create mode 100644 stabilizer.py diff --git a/stabilize.py b/stabilize.py index 29b0829..a40eef6 100644 --- a/stabilize.py +++ b/stabilize.py @@ -1,7 +1,8 @@ import SoapySDR import numpy as np from scipy.signal import blackmanharris -import serial + +from stabilizer import InductionHeater fs = 5e6 base_freq = 1086e6 @@ -9,84 +10,68 @@ threshold = 0.4 target_freq = 1088.3e6 k = 200e-6 -induction_min = 430e3 -induction_max = 445e3 - -generator = serial.Serial("/dev/ttyUSB0", 57600, writeTimeout=0) - -def set_induction_freq(freq): - command = ":s1f{:010d}\n".format(int(freq*1e2)) - generator.write(command.encode()) - generator.readline() - -def set_induction(amount): - amount = max(min(amount, 0.5), -0.5) - freq = (induction_min + induction_max)/2 + amount*(induction_max - induction_min) - set_induction_freq(freq) - - def parabolic(f, x): - """Quadratic interpolation for estimating the true position of an - inter-sample maximum when nearby samples are known. - f is a vector and x is an index for that vector. - Returns (vx, vy), the coordinates of the vertex of a parabola that goes - through point x and its two neighbors. - Example: - Defining a vector f with a local maximum at index 3 (= 6), find local - maximum if points 2, 3, and 4 actually defined a parabola. - In [3]: f = [2, 3, 1, 6, 4, 2, 3, 1] - In [4]: parabolic(f, argmax(f)) - Out[4]: (3.2142857142857144, 6.1607142857142856) - """ xv = 1/2. * (f[x-1] - f[x+1]) / (f[x-1] - 2 * f[x] + f[x+1]) + x yv = f[x] - 1/4. * (f[x-1] - f[x+1]) * (xv - x) return (xv, yv) + def fmt_range(x, m): pos = int(round(x*80/m)) return " "*pos + "*" -sdr = SoapySDR.Device() -samples = np.array([0]*4096, np.complex64) -sdr.setSampleRate(SoapySDR.SOAPY_SDR_RX, 0, fs) -sdr.setFrequency(SoapySDR.SOAPY_SDR_RX, 0, base_freq) +def main(): + induction = InductionHeater("/dev/ttyUSB0", 430e3, 445e3) + induction.start() + try: + sdr = SoapySDR.Device() + samples = np.array([0]*4096, np.complex64) -rxStream = sdr.setupStream(SoapySDR.SOAPY_SDR_RX, SoapySDR.SOAPY_SDR_CF32) -sdr.activateStream(rxStream) + sdr.setSampleRate(SoapySDR.SOAPY_SDR_RX, 0, fs) + sdr.setFrequency(SoapySDR.SOAPY_SDR_RX, 0, base_freq) -decimation = 0 + rxStream = sdr.setupStream(SoapySDR.SOAPY_SDR_RX, SoapySDR.SOAPY_SDR_CF32) + sdr.activateStream(rxStream) -try: - while True: - sr = sdr.readStream(rxStream, [samples], len(samples)) - if sr.ret != len(samples): - print("!") - continue + decimation = 0 - freqs = np.abs(np.fft.fft(samples*blackmanharris(len(samples)))[0:len(samples)//2]) - for i in range(len(freqs)//100): - freqs[i] = 0 - freqs[-i] = 0 + try: + while True: + sr = sdr.readStream(rxStream, [samples], len(samples)) + if sr.ret != len(samples): + print("!") + continue - # https://gist.github.com/endolith/255291 - i = np.argmax(freqs) - true_i = parabolic(freqs, i)[0] - freq = 0.5 * fs * true_i / len(freqs) - amplitude = freqs[i] + freqs = np.abs(np.fft.fft(samples*blackmanharris(len(samples)))[0:len(samples)//2]) + for i in range(len(freqs)//100): + freqs[i] = 0 + freqs[-i] = 0 - if amplitude > threshold: - actual_freq = base_freq + freq - print("{:8.3f} {:6.2f} {}".format(actual_freq/1e6, amplitude, fmt_range(freq, fs/2))) - induction = (actual_freq - target_freq)*k - else: - print("", amplitude) - induction = 0.0 + # https://gist.github.com/endolith/255291 + i = np.argmax(freqs) + true_i = parabolic(freqs, i)[0] + freq = 0.5 * fs * true_i / len(freqs) + amplitude = freqs[i] - decimation += 1 - if decimation == 15: - set_induction(induction) - decimation = 0 -finally: - sdr.deactivateStream(rxStream) - sdr.closeStream(rxStream) + if amplitude > threshold: + actual_freq = base_freq + freq + print("{:8.3f} {:6.2f} {}".format(actual_freq/1e6, amplitude, fmt_range(freq, fs/2))) + induction_amt = (actual_freq - target_freq)*k + else: + print("", amplitude) + induction_amt = 0.0 + + decimation += 1 + if decimation == 30: + induction.set(induction_amt) + decimation = 0 + finally: + sdr.deactivateStream(rxStream) + sdr.closeStream(rxStream) + finally: + induction.stop() + + +if __name__ == "__main__": + main() diff --git a/stabilizer.py b/stabilizer.py new file mode 100644 index 0000000..7fe2504 --- /dev/null +++ b/stabilizer.py @@ -0,0 +1,39 @@ +import serial +import queue +import threading + + +class InductionHeater: + """Interface to the MHS5200A function generator driving the LC tank""" + def __init__(self, port, induction_min, induction_max): + self.port = port + self.induction_min = induction_min + self.induction_max = induction_max + self.queue = queue.Queue(1) + + def start(self): + self.serial = serial.Serial(self.port, 57600) + self.thread = threading.Thread(target=self.thread_target) + self.thread.start() + + def thread_target(self): + while True: + amount = self.queue.get() + if amount is None: + break + + amount = max(min(amount, 0.5), -0.5) + freq = ((self.induction_min + self.induction_max)/2 + + amount*(self.induction_max - self.induction_min)) + + command = ":s1f{:010d}\n".format(int(freq*1e2)) + self.serial.write(command.encode()) + self.serial.readline() + + def set(self, amount): + self.queue.put(amount, block=False) + + def stop(self): + self.queue.put(None, block=True) + self.thread.join() + self.serial.close()