Mercurial > ~darius > hgwebdir.cgi > pyinst
comparison phasediff.py @ 61:b6dc43e6f6f0
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() |