complete DMI code, untested

dmi.py

@@ -0,0 +1,53 @@
+import SoapySDR
+import numpy as np
+
+from noptica import InductionHeater, Stabilizer, PositionTracker
+
+
+def main():
+    freq_sample = 5e6
+    freq_base = 1086e6
+    bufsize = 4096
+
+    induction = InductionHeater("/dev/ttyUSB0", 430e3, 445e3)
+    induction.start()
+    try:
+        stabilizer = Stabilizer(freq_sample, 0.4, 1088.3e6 - freq_base, 200e-6, induction)
+        position_tracker = PositionTracker(int(0.1*freq_sample/bufsize))
+
+        sdr = SoapySDR.Device()
+        for channel in range(2):
+            sdr.setSampleRate(SoapySDR.SOAPY_SDR_RX, channel, freq_sample)
+            sdr.setFrequency(SoapySDR.SOAPY_SDR_RX, channel, freq_base)
+
+        samples_ref = np.array([0]*bufsize, np.complex64)
+        samples_meas = np.array([0]*bufsize, np.complex64)
+        rx_stream = sdr.setupStream(SoapySDR.SOAPY_SDR_RX, SoapySDR.SOAPY_SDR_CF32, [0, 1])
+        try:
+            sdr.activateStream(rx_stream)
+            try:
+                stabilizer_throttle = 0
+                while True:
+                    sr = sdr.readStream(rx_stream, [samples_ref, samples_meas], bufsize)
+                    if sr.ret != bufsize:
+                        print("SDR sampling error")
+                        return
+
+                    # We can't update faster than the MHS5200A serial interface
+                    stabilizer_throttle += 1
+                    if stabilizer_throttle == 30:
+                        stabilizer.input(samples_ref)
+                        stabilizer_throttle = 0
+
+                    position, leakage = position_tracker.input(samples_ref, samples_meas)
+                    print(np.sum(position)/len(position), leakage)
+            finally:
+                sdr.deactivateStream(rx_stream)
+        finally:
+            sdr.closeStream(rx_stream)
+    finally:
+        induction.stop()
+
+
+if __name__ == "__main__":
+    main()

noptica.py

@@ -0,0 +1,102 @@
+import serial
+import queue
+import threading
+
+import numpy as np
+from scipy.signal import blackmanharris
+
+
+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()
+
+
+# https://gist.github.com/endolith/255291
+def parabolic(f, x):
+    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)
+
+
+class Stabilizer:
+    def __init__(self, freq_sample, amp_threshold, freq_target, k, tuner):
+        self.freq_sample = freq_sample
+        self.amp_threshold = amp_threshold
+        self.freq_target = freq_target
+        self.k = k
+        self.tuner = tuner
+
+    def input(self, samples):
+        spectrum = np.abs(np.fft.fft(samples*blackmanharris(len(samples)))[0:len(samples)//2])
+        for i in range(len(spectrum)//100):
+            spectrum[i] = 0
+            spectrum[-i] = 0
+        i = np.argmax(spectrum)
+        true_i, amplitude = parabolic(spectrum, i)
+        freq = 0.5 * self.freq_sample * true_i / len(spectrum)
+
+        if amplitude > threshold:
+            tuning = (freq - target_freq)*k
+        else:
+            tuning = 0.0
+        self.tuner.set(tuning)
+
+
+def continuous_unwrap(last_phase, last_phase_unwrapped, p):
+    # note: np.unwrap always preserves first element of array
+    p = np.unwrap(p)
+    glue = np.array([last_phase_unwrapped, last_phase_unwrapped + (p[0] - last_phase)])
+    new_p0 = np.unwrap(glue)[1]
+    return new_p0 + p - p[0]
+
+
+class PositionTracker:
+    def __init__(self, leakage_avg):
+        self.last_phase = 0.0
+        self.last_position = 0.0
+        self.leakage = np.zeros(leakage_avg)
+        self.leakage_ptr = 0
+
+    def input(self, ref, meas):
+        demod = np.conjugate(ref)*meas
+
+        self.leakage[self.leakage_ptr] = np.real(np.sum(demod)/len(demod))
+        self.leakage_ptr = (self.leakage_ptr + 1) % len(self.leakage)
+        leakage = np.sum(self.leakage)/len(self.leakage)
+
+        phase = np.angle(demod - leakage)
+        position = continuous_unwrap(self.last_phase, self.last_position, phase)/(2.0*np.pi)
+        self.last_phase = phase[-1]
+        self.last_position = position[-1]
+        return position, leakage

stabilize.py

@@ -1,77 +0,0 @@
-import SoapySDR
-import numpy as np
-from scipy.signal import blackmanharris
-
-from stabilizer import InductionHeater
-
-fs = 5e6
-base_freq = 1086e6
-threshold = 0.4
-target_freq = 1088.3e6
-k = 200e-6
-
-def parabolic(f, x):
-    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 + "*"
-
-
-def main():
-    induction = InductionHeater("/dev/ttyUSB0", 430e3, 445e3)
-    induction.start()
-    try:
-        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)
-
-        rxStream = sdr.setupStream(SoapySDR.SOAPY_SDR_RX, SoapySDR.SOAPY_SDR_CF32)
-        sdr.activateStream(rxStream)
-
-        decimation = 0
-
-        try:
-            while True:
-                sr = sdr.readStream(rxStream, [samples], len(samples))
-                if sr.ret != len(samples):
-                    print("!")
-                    continue
-
-                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
-
-                # 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]
-
-                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("<no signal>", 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()

stabilizer.py

@@ -1,39 +0,0 @@
-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()