# HG changeset patch # User Daniel O'Connor # Date 1610076987 -37800 # Node ID b6dc43e6f6f06773eadb3ce8cb9e44228c6bdf91 # Parent 4558e5ccd775127948da17af6fa71002dc6d197c Add code to measure phase difference between 2 signals by sampling from the cro. diff -r 4558e5ccd775 -r b6dc43e6f6f0 phasediff.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/phasediff.py Fri Jan 08 14:06:27 2021 +1030 @@ -0,0 +1,230 @@ +#!/usr/bin/env python + + + +# Copyright (c) 2020 +# Daniel O'Connor . All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# 1. Redistributions of source code must retain the above copyright +# notice, this list of conditions and the following disclaimer. +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND +# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE +# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +# OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY +# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF +# SUCH DAMAGE. +# + +# Expected DB schema +# CREATE TABLE phasediff ( +# name TEXT, +# time TIMESTAMP WITH TIME ZONE, +# delta NUMERIC(15, 12), +# cumulative NUMERIC(15, 12) +# ); + +import datetime +import matplotlib.pylab as pylab +import numpy +import psycopg2 +import scipy.signal +import sqlite3 +import sys +import time +# Should try this code instead: https://github.com/python-ivi/python-usbtmc +import usb488 + +def main(): + u = usb488.USB488Device() + print('Found device') + + # See "TDS2000 Programmer.pdf" + u.write('*IDN?') + print(('IDN reports ' + u.read())) + + #dbh = sqlite3.connect('phasediff.db') + dbh = psycopg2.connect('host=vm11 user=phasediff dbname=phasediff') + test(u, dbh, 'delamerelts') + +def test(u, dbh = None, name = None): + if dbh != None: + cur = dbh.cursor() + + u.write('DATA:ENC RIB') # Big endian signed + u.write('DATA:WIDTH 2') # 2 bytes wide + + u.write('CH1:SCALE?') + vscale1 = float(u.read(1).split()[1]) + print(('Channel 1 scale is %.2f volts/div' % (vscale1))) + u.write('CH2:SCALE?') + vscale2 = float(u.read(1).split()[1]) + print(('Channel 2 scale is %.2f volts/div' % (vscale2))) + u.write('HOR:MAIN:SCALE?') + hscale = float(u.read(1).split()[1]) + print(('Horizontal scale is %.5f nsec/div' % (hscale * 1e9))) + # TEK2024B doesn't grok HOR:DIV? so hard code 10 (has 8 vertically) + acqwindow = hscale * 10.0 + nomclockperiod = 1 / 10.7e6 + + cumdiff = 0 + lastdiff = None + while True: + ary1, ary2 = acquire(u, vscale1, vscale2) + + sampletime = acqwindow / len(ary1) + + # Compute clock periods + # XXX: this is pretty noisy so we use the nominal clock period above for unwrapping + #period1 = getperiod(ary1) * sampletime + #period2 = getperiod(ary2) * sampletime + + # Compute phased difference between the two + d = getpdiff(ary1, ary2) * sampletime + #d = getpdiffedge(ary1, ary2) * sampletime + if lastdiff == None: + # First time, init variables and start again + # XXX: should get the most recent cumuulative diff from the DB really.. + lastdiff = d + cumdiff = d + continue + + # Difference between difference + diff = d - lastdiff + + # Check if it wrapped + if abs(diff) > nomclockperiod / 2: + if d > nomclockperiod / 2: + diff -= nomclockperiod + else: + diff += nomclockperiod + + # Accumulate deltas + cumdiff += diff + + # Log / insert into DB + now = datetime.datetime.now() + print(("%s Cumulative: %.2f nsec Delta: %.2f nsec Diff: %.2f nsec" % ( + now.strftime('%Y/%m/%d %H:%M:%S.%f'), + cumdiff * 1e9, d * 1e9, diff * 1e9))) + if dbh != None: + cur.execute('INSERT INTO phasediff(name, time, delta, cumulative) VALUES(%s, %s, %s, %s)', (name, now, d, cumdiff)) + dbh.commit() + lastdiff = d + +def acquire(u, vscale1, vscale2): + u.write('ACQ:STATE 1') # Do a single acquisition + u.write('*OPC?') + u.read(2.0) # Wait for it to complete + + u.write('DAT:SOU CH1') # Set the curve source to channel 1 + u.write('CURVE?') # Ask for the curve data + then = time.time() + result = u.read(1.0) # Takes the CRO a while for this + now = time.time() + #print('CURVE read took %f milliseconds' % ((now - then) * 1000.0)) + data1 = result[13:] # Chop off the header (should verify this really..) + ary1 = numpy.fromstring(data1, dtype = '>h') + ary1 = ary1 / 32768.0 * vscale1 # Scale to volts + + u.write('DAT:SOU CH2') # Set the curve source to channel 2 + u.write('CURVE?') # Ask for the curve data + then = time.time() + result = u.read(1.0) # Takes the CRO a while for this + now = time.time() + #print('CURVE read took %f milliseconds' % ((now - then) * 1000.0)) + data2 = result[13:] # Chop off the header (should verify this really..) + ary2 = numpy.fromstring(data2, dtype = '>h') + ary2 = ary2 / 32768.0 * vscale2 # Scale to volts + + return ary1, ary2 + +def getpdiff(ary1, ary2): + '''Return phase difference in samples between two signals by cross correlation''' + + # Rescale to 0-1 + ary1 = ary1 - ary1.min() + ary1 = ary1 / ary1.max() + ary2 = ary2 - ary2.min() + ary2 = ary2 / ary2.max() + + # Cross correlate + corr = scipy.signal.correlate(ary1, ary2) + + # Find peak + amax = corr.argmax() + + return amax - (len(ary1) - 1) + +def getpdiffedge(ary1, ary2): + '''Return phase difference in samples between two signals by edge detection''' + + # Rescale to 0-1 + ary1 = ary1 - ary1.min() + ary1 = ary1 / ary1.max() + ary2 = ary2 - ary2.min() + ary2 = ary2 / ary2.max() + + # Find rising edge of each + ary1pos = numpy.argmax(ary1 > 0.2) + ary2pos = numpy.argmax(ary2 > 0.2) + + return ary1pos - ary2pos + +def getperiod(ary): + '''Return period of signal in ary in samples''' + # Compute auto correlation + corr = scipy.signal.correlate(ary, ary) + + # Find peaks + peaks, _ = scipy.signal.find_peaks(corr) + + # Find time differences between peaks + deltapeak = (peaks[1:-1] - peaks[0:-2]) + + # Return average of the middle few + return deltapeak[2:-2].mean() + +def trend(dbh, name): + cur = dbh.cursor() + cur2 = dbh.cursor() + cur.execute('SELECT time, deltansec FROM phasediff WHERE name = %s', (name, )) + + _, lastdelta = cur.fetchone() + + count = 0 + cumdiff = 0 + for row in cur: + time, delta = row + diff = float(delta - lastdelta) + if abs(diff) > 45e-9: + if delta > 45e-9: + diff -= 1.0 / 10.7e6 + else: + diff += 1.0 / 10.7e6 + cumdiff += diff + #print('Cumulative: %.5f nsec Delta: %.5f nsec Last: %.5f nsec Diff: %.5f nsec' % ( + # float(cumdiff) * 1e9, float(delta) * 1e9, float(lastdelta) * 1e9, float(diff) * 1e9)) + lastdelta = delta + cur2.execute('INSERT INTO cumdiff (name, time, cummdiff) VALUES (%s, %s, %s)', + ( name, time, cumdiff)) + count += 1 + #if count % 20 == 0: + # sys.stdin.readline() + if count % 1000 == 0: + dbh.commit() + print(('Count %d: %s' % (count, str(time)))) + +if __name__ == '__main__': + main()