Mercurial > ~darius > hgwebdir.cgi > pyinst

comparison phasediff.py @ 61:b6dc43e6f6f0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add code to measure phase difference between 2 signals by sampling from the cro.
author Daniel O'Connor <doconnor@gsoft.com.au>
date Fri, 08 Jan 2021 14:06:27 +1030
parents
children
comparison
equal deleted inserted replaced
60:4558e5ccd775 61:b6dc43e6f6f0
1 #!/usr/bin/env python
2
3
4
5 # Copyright (c) 2020
6 # Daniel O'Connor <darius@dons.net.au>. All rights reserved.
7 #
8 # Redistribution and use in source and binary forms, with or without
9 # modification, are permitted provided that the following conditions
10 # are met:
11 # 1. Redistributions of source code must retain the above copyright
12 # notice, this list of conditions and the following disclaimer.
13 # 2. Redistributions in binary form must reproduce the above copyright
14 # notice, this list of conditions and the following disclaimer in the
15 # documentation and/or other materials provided with the distribution.
16 #
17 # THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND
18 # ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 # ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
21 # FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
22 # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
23 # OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
24 # HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
25 # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
26 # OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
27 # SUCH DAMAGE.
28 #
29
30 # Expected DB schema
31 # CREATE TABLE phasediff (
32 # name TEXT,
33 # time TIMESTAMP WITH TIME ZONE,
34 # delta NUMERIC(15, 12),
35 # cumulative NUMERIC(15, 12)
36 # );
37
38 import datetime
39 import matplotlib.pylab as pylab
40 import numpy
41 import psycopg2
42 import scipy.signal
43 import sqlite3
44 import sys
45 import time
46 # Should try this code instead: https://github.com/python-ivi/python-usbtmc
47 import usb488
48
49 def main():
50 u = usb488.USB488Device()
51 print('Found device')
52
53 # See "TDS2000 Programmer.pdf"
54 u.write('*IDN?')
55 print(('IDN reports ' + u.read()))
56
57 #dbh = sqlite3.connect('phasediff.db')
58 dbh = psycopg2.connect('host=vm11 user=phasediff dbname=phasediff')
59 test(u, dbh, 'delamerelts')
60
61 def test(u, dbh = None, name = None):
62 if dbh != None:
63 cur = dbh.cursor()
64
65 u.write('DATA:ENC RIB') # Big endian signed
66 u.write('DATA:WIDTH 2') # 2 bytes wide
67
68 u.write('CH1:SCALE?')
69 vscale1 = float(u.read(1).split()[1])
70 print(('Channel 1 scale is %.2f volts/div' % (vscale1)))
71 u.write('CH2:SCALE?')
72 vscale2 = float(u.read(1).split()[1])
73 print(('Channel 2 scale is %.2f volts/div' % (vscale2)))
74 u.write('HOR:MAIN:SCALE?')
75 hscale = float(u.read(1).split()[1])
76 print(('Horizontal scale is %.5f nsec/div' % (hscale * 1e9)))
77 # TEK2024B doesn't grok HOR:DIV? so hard code 10 (has 8 vertically)
78 acqwindow = hscale * 10.0
79 nomclockperiod = 1 / 10.7e6
80
81 cumdiff = 0
82 lastdiff = None
83 while True:
84 ary1, ary2 = acquire(u, vscale1, vscale2)
85
86 sampletime = acqwindow / len(ary1)
87
88 # Compute clock periods
89 # XXX: this is pretty noisy so we use the nominal clock period above for unwrapping
90 #period1 = getperiod(ary1) * sampletime
91 #period2 = getperiod(ary2) * sampletime
92
93 # Compute phased difference between the two
94 d = getpdiff(ary1, ary2) * sampletime
95 #d = getpdiffedge(ary1, ary2) * sampletime
96 if lastdiff == None:
97 # First time, init variables and start again
98 # XXX: should get the most recent cumuulative diff from the DB really..
99 lastdiff = d
100 cumdiff = d
101 continue
102
103 # Difference between difference
104 diff = d - lastdiff
105
106 # Check if it wrapped
107 if abs(diff) > nomclockperiod / 2:
108 if d > nomclockperiod / 2:
109 diff -= nomclockperiod
110 else:
111 diff += nomclockperiod
112
113 # Accumulate deltas
114 cumdiff += diff
115
116 # Log / insert into DB
117 now = datetime.datetime.now()
118 print(("%s Cumulative: %.2f nsec Delta: %.2f nsec Diff: %.2f nsec" % (
119 now.strftime('%Y/%m/%d %H:%M:%S.%f'),
120 cumdiff * 1e9, d * 1e9, diff * 1e9)))
121 if dbh != None:
122 cur.execute('INSERT INTO phasediff(name, time, delta, cumulative) VALUES(%s, %s, %s, %s)', (name, now, d, cumdiff))
123 dbh.commit()
124 lastdiff = d
125
126 def acquire(u, vscale1, vscale2):
127 u.write('ACQ:STATE 1') # Do a single acquisition
128 u.write('*OPC?')
129 u.read(2.0) # Wait for it to complete
130
131 u.write('DAT:SOU CH1') # Set the curve source to channel 1
132 u.write('CURVE?') # Ask for the curve data
133 then = time.time()
134 result = u.read(1.0) # Takes the CRO a while for this
135 now = time.time()
136 #print('CURVE read took %f milliseconds' % ((now - then) * 1000.0))
137 data1 = result[13:] # Chop off the header (should verify this really..)
138 ary1 = numpy.fromstring(data1, dtype = '>h')
139 ary1 = ary1 / 32768.0 * vscale1 # Scale to volts
140
141 u.write('DAT:SOU CH2') # Set the curve source to channel 2
142 u.write('CURVE?') # Ask for the curve data
143 then = time.time()
144 result = u.read(1.0) # Takes the CRO a while for this
145 now = time.time()
146 #print('CURVE read took %f milliseconds' % ((now - then) * 1000.0))
147 data2 = result[13:] # Chop off the header (should verify this really..)
148 ary2 = numpy.fromstring(data2, dtype = '>h')
149 ary2 = ary2 / 32768.0 * vscale2 # Scale to volts
150
151 return ary1, ary2
152
153 def getpdiff(ary1, ary2):
154 '''Return phase difference in samples between two signals by cross correlation'''
155
156 # Rescale to 0-1
157 ary1 = ary1 - ary1.min()
158 ary1 = ary1 / ary1.max()
159 ary2 = ary2 - ary2.min()
160 ary2 = ary2 / ary2.max()
161
162 # Cross correlate
163 corr = scipy.signal.correlate(ary1, ary2)
164
165 # Find peak
166 amax = corr.argmax()
167
168 return amax - (len(ary1) - 1)
169
170 def getpdiffedge(ary1, ary2):
171 '''Return phase difference in samples between two signals by edge detection'''
172
173 # Rescale to 0-1
174 ary1 = ary1 - ary1.min()
175 ary1 = ary1 / ary1.max()
176 ary2 = ary2 - ary2.min()
177 ary2 = ary2 / ary2.max()
178
179 # Find rising edge of each
180 ary1pos = numpy.argmax(ary1 > 0.2)
181 ary2pos = numpy.argmax(ary2 > 0.2)
182
183 return ary1pos - ary2pos
184
185 def getperiod(ary):
186 '''Return period of signal in ary in samples'''
187 # Compute auto correlation
188 corr = scipy.signal.correlate(ary, ary)
189
190 # Find peaks
191 peaks, _ = scipy.signal.find_peaks(corr)
192
193 # Find time differences between peaks
194 deltapeak = (peaks[1:-1] - peaks[0:-2])
195
196 # Return average of the middle few
197 return deltapeak[2:-2].mean()
198
199 def trend(dbh, name):
200 cur = dbh.cursor()
201 cur2 = dbh.cursor()
202 cur.execute('SELECT time, deltansec FROM phasediff WHERE name = %s', (name, ))
203
204 _, lastdelta = cur.fetchone()
205
206 count = 0
207 cumdiff = 0
208 for row in cur:
209 time, delta = row
210 diff = float(delta - lastdelta)
211 if abs(diff) > 45e-9:
212 if delta > 45e-9:
213 diff -= 1.0 / 10.7e6
214 else:
215 diff += 1.0 / 10.7e6
216 cumdiff += diff
217 #print('Cumulative: %.5f nsec Delta: %.5f nsec Last: %.5f nsec Diff: %.5f nsec' % (
218 # float(cumdiff) * 1e9, float(delta) * 1e9, float(lastdelta) * 1e9, float(diff) * 1e9))
219 lastdelta = delta
220 cur2.execute('INSERT INTO cumdiff (name, time, cummdiff) VALUES (%s, %s, %s)',
221 ( name, time, cumdiff))
222 count += 1
223 #if count % 20 == 0:
224 # sys.stdin.readline()
225 if count % 1000 == 0:
226 dbh.commit()
227 print(('Count %d: %s' % (count, str(time))))
228
229 if __name__ == '__main__':
230 main()