Merge branch 'wip' of gitlab.psi.ch-samenv:samenv/frappy into wip

This commit is contained in:
zolliker 2022-05-30 12:10:51 +02:00
commit 4cb0031302
4 changed files with 437 additions and 32 deletions

294
secop_psi/adq_mr.py Normal file
View File

@ -0,0 +1,294 @@
# -*- coding: utf-8 -*-
"""
Created on Tue Nov 26 15:42:43 2019
@author: tartarotti_d-adm
"""
import numpy as np
import ctypes as ct
import time
from numpy import sqrt, arctan2, sin, cos
#from pylab import *
from scipy import signal
#ADQAPI = ct.cdll.LoadLibrary("ADQAPI.dll")
ADQAPI = ct.cdll.LoadLibrary("libadq.so.0")
#For different trigger modes
SW_TRIG = 1
EXT_TRIG_1 = 2 #This external trigger does not work if the level of the trigger is very close to 0.5V. Now we have it close to 3V, and it works
EXT_TRIG_2 = 7
EXT_TRIG_3 = 8
LVL_TRIG = 3
INT_TRIG = 4
LVL_FALLING = 0
LVL_RISING = 1
#samples_per_record=16384
ADQ_TRANSFER_MODE_NORMAL = 0x00
ADQ_CHANNELS_MASK = 0x3
#f_LO = 40
def butter_lowpass(cutoff, sr, order=5):
nyq = 0.5 * sr
normal_cutoff = cutoff / nyq
b, a = signal.butter(order, normal_cutoff, btype = 'low', analog = False)
return b, a
class Adq(object):
max_number_of_channels = 2
samp_freq = 2
#ndecimate = 50 # decimation ratio (2GHz / 40 MHz)
ndecimate = 50
def __init__(self, number_of_records, samples_per_record, bw_cutoff):
self.number_of_records = number_of_records
self.samples_per_record = samples_per_record
self.bw_cutoff = bw_cutoff
ADQAPI.ADQAPI_GetRevision()
# Manually set return type from some ADQAPI functions
ADQAPI.CreateADQControlUnit.restype = ct.c_void_p
ADQAPI.ADQ_GetRevision.restype = ct.c_void_p
ADQAPI.ADQ_GetPtrStream.restype = ct.POINTER(ct.c_int16)
ADQAPI.ADQControlUnit_FindDevices.argtypes = [ct.c_void_p]
# Create ADQControlUnit
self.adq_cu = ct.c_void_p(ADQAPI.CreateADQControlUnit())
ADQAPI.ADQControlUnit_EnableErrorTrace(self.adq_cu, 3, '.')
self.adq_num = 1
# Find ADQ devices
ADQAPI.ADQControlUnit_FindDevices(self.adq_cu)
n_of_ADQ = ADQAPI.ADQControlUnit_NofADQ(self.adq_cu)
if n_of_ADQ != 1:
raise ValueError('number of ADQs must be 1, not %d' % n_of_ADQ)
rev = ADQAPI.ADQ_GetRevision(self.adq_cu, self.adq_num)
revision = ct.cast(rev,ct.POINTER(ct.c_int))
print('\nConnected to ADQ #1')
# Print revision information
print('FPGA Revision: {}'.format(revision[0]))
if (revision[1]):
print('Local copy')
else :
print('SVN Managed')
if (revision[2]):
print('Mixed Revision')
else :
print('SVN Updated')
print('')
ADQ_CLOCK_INT_INTREF = 0 #internal clock source
ADQ_CLOCK_EXT_REF = 1 #internal clock source, external reference
ADQ_CLOCK_EXT_CLOCK = 2 #External clock source
ADQAPI.ADQ_SetClockSource(self.adq_cu, self.adq_num, ADQ_CLOCK_EXT_REF);
##########################
# Test pattern
#ADQAPI.ADQ_SetTestPatternMode(self.adq_cu, self.adq_num, 4)
##########################
# Sample skip
#ADQAPI.ADQ_SetSampleSkip(self.adq_cu, self.adq_num, 1)
##########################
# Set trig mode
self.trigger = EXT_TRIG_1
#trigger = LVL_TRIG
success = ADQAPI.ADQ_SetTriggerMode(self.adq_cu, self.adq_num, self.trigger)
if (success == 0):
print('ADQ_SetTriggerMode failed.')
if (self.trigger == LVL_TRIG):
success = ADQAPI.ADQ_SetLvlTrigLevel(self.adq_cu, self.adq_num, -100)
if (success == 0):
print('ADQ_SetLvlTrigLevel failed.')
success = ADQAPI.ADQ_SetTrigLevelResetValue(self.adq_cu, self.adq_num, 1000)
if (success == 0):
print('ADQ_SetTrigLevelResetValue failed.')
success = ADQAPI.ADQ_SetLvlTrigChannel(self.adq_cu, self.adq_num, 1)
if (success == 0):
print('ADQ_SetLvlTrigChannel failed.')
success = ADQAPI.ADQ_SetLvlTrigEdge(self.adq_cu, self.adq_num, LVL_RISING)
if (success == 0):
print('ADQ_SetLvlTrigEdge failed.')
elif (self.trigger == EXT_TRIG_1) :
success = ADQAPI.ADQ_SetExternTrigEdge(self.adq_cu, self.adq_num,2)
if (success == 0):
print('ADQ_SetLvlTrigEdge failed.')
# success = ADQAPI.ADQ_SetTriggerThresholdVoltage(self.adq_cu, self.adq_num, trigger, ct.c_double(0.2))
# if (success == 0):
# print('SetTriggerThresholdVoltage failed.')
print("CHANNEL:"+str(ct.c_int(ADQAPI.ADQ_GetLvlTrigChannel(self.adq_cu, self.adq_num))))
self.setup_target_buffers()
def setup_target_buffers(self):
# Setup target buffers for data
self.target_buffers=(ct.POINTER(ct.c_int16 * self.samples_per_record * self.number_of_records)
* self.max_number_of_channels)()
for bufp in self.target_buffers:
bufp.contents = (ct.c_int16 * self.samples_per_record * self.number_of_records)()
def deletecu(self):
# Only disarm trigger after data is collected
ADQAPI.ADQ_DisarmTrigger(self.adq_cu, self.adq_num)
ADQAPI.ADQ_MultiRecordClose(self.adq_cu, self.adq_num);
# Delete ADQControlunit
ADQAPI.DeleteADQControlUnit(self.adq_cu)
def start(self):
"""start datat acquisition"""
# samples_per_records = samples_per_record/number_of_records
# Change number of pulses to be acquired acording to how many records are taken
# Start acquisition
ADQAPI.ADQ_MultiRecordSetup(self.adq_cu, self.adq_num,
self.number_of_records,
self.samples_per_record)
ADQAPI.ADQ_DisarmTrigger(self.adq_cu, self.adq_num)
ADQAPI.ADQ_ArmTrigger(self.adq_cu, self.adq_num)
def getdata(self):
"""wait for aquisition to be finished and get data"""
#start = time.time()
while(ADQAPI.ADQ_GetAcquiredAll(self.adq_cu,self.adq_num) == 0):
time.sleep(0.001)
#if (self.trigger == SW_TRIG):
# ADQAPI.ADQ_SWTrig(self.adq_cu, self.adq_num)
#mid = time.time()
status = ADQAPI.ADQ_GetData(self.adq_cu, self.adq_num, self.target_buffers,
self.samples_per_record * self.number_of_records, 2,
0, self.number_of_records, ADQ_CHANNELS_MASK,
0, self.samples_per_record, ADQ_TRANSFER_MODE_NORMAL);
#print(time.time()-mid,mid-start)
if not status:
raise ValueError('no succesS from ADQ_GetDATA')
# Now this is an array with all records, but the time is artificial
data = []
for ch in range(2):
onedim = np.frombuffer(self.target_buffers[ch].contents, dtype=np.int16)
data.append(onedim.reshape(self.number_of_records, self.samples_per_record) / float(2**14)) # 14 bits ADC
return data
def acquire(self):
self.start()
return self.getdata()
'''
def average(self, data):
#Average over records
return [data[ch].sum(axis=0) / self.number_of_records for ch in range(2)]
def iq(self, channel, f_LO):
newx = np.linspace(0, self.samples_per_record /2, self.samples_per_record)
s0 = channel /((2**16)/2)*0.5*np.exp(1j*2*np.pi*f_LO/(1e3)*newx)
I0 = s0.real
Q0 = s0.imag
return I0, Q0
def fitting(self, data, f_LO, ti, tf):
# As long as data[0] is the pulse
si = 2*ti #Those are for fitting the pulse
sf = 2*tf
phase = np.zeros(self.number_of_records)
amplitude = np.zeros(self.number_of_records)
offset = np.zeros(self.number_of_records)
for i in range(self.number_of_records):
phase[i], amplitude[i] = sineW(data[0][i][si:sf],f_LO*1e-9,ti,tf)
offset[i] = np.average(data[0][i][si:sf])
return phase, amplitude, offset
def waveIQ(self, channel,ti,f_LO):
#channel is not the sample data
t = np.linspace(0, self.samples_per_record /2, self.samples_per_record + 1)[:-1]
si = 2*ti # Again that is where the wave pulse starts
cwi = np.zeros((self.number_of_records,self.samples_per_record))
cwq = np.zeros((self.number_of_records,self.samples_per_record))
iq = np.zeros((self.number_of_records,self.samples_per_record))
q = np.zeros((self.number_of_records,self.samples_per_record))
for i in range(self.number_of_records):
cwi[i] = np.zeros(self.samples_per_record)
cwq[i] = np.zeros(self.samples_per_record)
cwi[i] = amplitude[i]*sin(t[si:]*f_LO*1e-9*2*np.pi+phase[i]*np.pi/180)+bias[i]
cwq[i] = amplitude[i]*sin(t[si:]*f_LO*1e-9*(2*np.pi+(phase[i]+90)*np.pi/180))+bias[i]
iq[i] = channel[i]*cwi[i]
q[i] = channel[i]*cwq[i]
return iq,q
'''
def sinW(self,sig,freq,ti,tf):
# sig: signal array
# freq
# ti, tf: initial and end time
si = int(ti * self.samp_freq)
nperiods = freq * (tf - ti)
n = int(round(max(2, int(nperiods)) / nperiods * (tf-ti) * self.samp_freq))
self.nperiods = n
t = np.arange(si, len(sig)) / self.samp_freq
t = t[:n]
self.pulselen = n / self.samp_freq
sig = sig[si:si+n]
a = 2*np.sum(sig*np.cos(2*np.pi*freq*t))/len(sig)
b = 2*np.sum(sig*np.sin(2*np.pi*freq*t))/len(sig)
return a, b
def mix(self, sigin, sigout, freq, ti, tf):
# sigin, sigout: signal array, incomping, output
# freq
# ti, tf: initial and end time if sigin
a, b = self.sinW(sigin, freq, ti, tf)
phase = arctan2(a,b) * 180 / np.pi
amp = sqrt(a**2 + b**2)
a, b = a/amp, b/amp
#si = int(ti * self.samp_freq)
t = np.arange(len(sigout)) / self.samp_freq
wave1 = sigout * (a * cos(2*np.pi*freq*t) + b * sin(2*np.pi*freq*t))
wave2 = sigout * (a * sin(2*np.pi*freq*t) - b * cos(2*np.pi*freq*t))
return wave1, wave2
def averageiq(self, data, freq, ti, tf):
'''Average over records'''
iorq = np.array([self.mix(data[0][i], data[1][i], freq, ti, tf) for i in range(self.number_of_records)])
# iorq = np.array([self.mix(data[0][:], data[1][:], freq, ti, tf)])
return iorq.sum(axis=0) / self.number_of_records
def filtro(self, iorq, cutoff):
b, a = butter_lowpass(cutoff, self.samp_freq*1e9)
#ifi = np.array(signal.filtfilt(b,a,iorq[0]))
#qf = np.array(signal.filtfilt(b,a,iorq[1]))
iqf = [signal.filtfilt(b,a,iorq[i]) for i in np.arange(len(iorq))]
return iqf
def box(self, iorq, ti, tf):
si = int(self.samp_freq * ti)
sf = int(self.samp_freq * tf)
bxa = [sum(iorq[i][si:sf])/(sf-si) for i in np.arange(len(iorq))]
return bxa
def gates_and_curves(self, data, freq, pulse, roi):
"""return iq values of rois and prepare plottable curves for iq"""
times = []
times.append(('aviq', time.time()))
iq = self.averageiq(data,freq*1e-9,*pulse)
times.append(('filtro', time.time()))
iqf = self.filtro(iq,self.bw_cutoff)
m = len(iqf[0]) // self.ndecimate
times.append(('iqdec', time.time()))
iqd = np.average(np.resize(iqf, (2, m, self.ndecimate)), axis=2)
t_axis = np.arange(m) * self.ndecimate / self.samp_freq
pulsig = np.abs(data[0][0])
times.append(('pulsig', time.time()))
pulsig = np.average(np.resize(pulsig, (m, self.ndecimate)), axis=1)
self.curves = (t_axis, iqd[0], iqd[1], pulsig)
#print(times)
return [self.box(iqf,*r) for r in roi]

View File

@ -20,43 +20,43 @@
# *****************************************************************************
"""Delay generator stanford 645"""
from secop.core import FloatRange, HasIodev, Module, Parameter, StringIO
from secop.core import FloatRange, HasIO, Module, Parameter, StringIO
class DG645(StringIO):
end_of_line = '\n'
class Delay(HasIodev, Module):
class Delay(HasIO, Module):
on1 = Parameter('on delay 1', FloatRange(unit='sec'), readonly=False, default=0)
off1 = Parameter('off delay 1', FloatRange(unit='sec'), readonly=False, default=60e-9)
on2 = Parameter('on delay 2', FloatRange(unit='sec'), readonly=False, default=0)
off2 = Parameter('off delay 2', FloatRange(unit='sec'), readonly=False, default=150e-9)
iodevClass = DG645
ioClass = DG645
def read_on1(self):
return self.sendRecv('DLAY?2').split(',')[1]
return float(self.communicate('DLAY?2').split(',')[1])
def read_off1(self):
return self.sendRecv('DLAY?3').split(',')[1]
return float(self.communicate('DLAY?3').split(',')[1])
def read_on2(self):
return self.sendRecv('DLAY?4').split(',')[1]
return float(self.communicate('DLAY?4').split(',')[1])
def read_off2(self):
return self.sendRecv('DLAY?5').split(',')[1]
return float(self.communicate('DLAY?5').split(',')[1])
def write_on1(self, value):
return self.sendRecv('DLAY 2,0,%g;DLAY?2' % value).split(',')[1]
return float(self.communicate('DLAY 2,0,%g;DLAY?2' % value).split(',')[1])
def write_off1(self, value):
result = self.sendRecv('DLAY 3,0,%g;DLAY?3' % value)
return result.split(',')[1]
result = self.communicate('DLAY 3,0,%g;DLAY?3' % value)
return float(result.split(',')[1])
def write_on2(self, value):
return self.sendRecv('DLAY 4,0,%g;DLAY?4' % value).split(',')[1]
return float(self.communicate('DLAY 4,0,%g;DLAY?4' % value).split(',')[1])
def write_off2(self, value):
return self.sendRecv('DLAY 5,0,%g;DLAY?5' % value).split(',')[1]
return float(self.communicate('DLAY 5,0,%g;DLAY?5' % value).split(',')[1])

102
secop_psi/iqplot.py Normal file
View File

@ -0,0 +1,102 @@
# -*- coding: utf-8 -*-
"""
Created on Tue Feb 4 11:07:56 2020
@author: tartarotti_d-adm
"""
import numpy as np
import matplotlib.pyplot as plt
def rect(x1, x2, y1, y2):
return np.array([[x1,x2,x2,x1,x1],[y1,y1,y2,y2,y1]])
NAN = float('nan')
def rects(intervals, y12):
result = [rect(*intervals[0], *y12)]
for x12 in intervals[1:]:
result.append([[NAN],[NAN]])
result.append(rect(*x12, *y12))
return np.concatenate(result, axis=1)
class Plot:
def __init__(self, maxy):
self.lines = {}
self.yaxis = ((-2 * maxy, maxy), (-maxy, 2 * maxy))
self.first = True
self.fig = None
def set_line(self, iax, name, data, fmt, **kwds):
"""
plot or update a line
when called with self.first = True: plot the line
when called with self.first = False: update the line
iax: 0: left, 1: right yaxis
name: the name of the line. used also as label for legend, if not starting with underscore
data: data[0] and data[1] are used for x/y data respectively
fmt, other keywords: forwarded to <axis>.plot
"""
# ax: 0: left, 1: right
if self.first:
if name.startswith('_'):
label = '_nolegend_'
else:
label = name
self.lines[name], = self.ax[iax].plot(data[0], data[1], fmt, label=label, **kwds)
else:
self.lines[name].set_data(data[0:2])
def close(self):
if self.fig:
plt.close(self.fig)
self.fig = None
self.first = True
def plot(self, curves, rois=None, average=None):
boxes = rects(rois[1:], self.yaxis[0])
pbox = rect(*rois[0], *self.yaxis[1])
rbox = rect(*rois[1], *self.yaxis[0])
pshift = self.yaxis[0][1] * 0.5
pulse = curves[3] - pshift
# normalized to 0.8 * pshift:
#pulse = (curves[3] / np.max(curves[3]))* pshift * 0.8 - pshift
try:
if self.first:
plt.ion()
self.fig, axleft = plt.subplots(figsize=(15,7))
plt.title("I/Q", fontsize=14)
axleft.set_xlim(0, curves[0][-1])
self.ax = [axleft, axleft.twinx()]
self.ax[0].axhline(y=0, color='#cccccc') # show x-axis line
self.ax[1].axhline(y=0, color='#cccccc')
self.ax[0].set_ylim(*self.yaxis[0])
self.ax[1].set_ylim(*self.yaxis[1])
self.set_line(0, "I", curves, 'b-') # using curves [0] and [1]
self.set_line(0, "_Iaverage", average, 'b.')
self.set_line(0, "Ampl", (curves[0],np.sqrt(curves[1]**2+curves[2]**2)), '#808080')
self.set_line(1, "Q", (curves[0], curves[2]), 'g-')
self.set_line(1, "_Qaverage", (average[0], average[2]), 'g.')
self.set_line(0, "pulse", (curves[0], pulse), 'c-')
self.set_line(0, "roi's", boxes, 'm-')
self.set_line(1, "pulse reg", pbox, 'k-')
self.set_line(0, "ctrl reg", rbox, 'r-')
if self.first:
self.fig.legend(fontsize=12)
plt.tight_layout()
finally:
self.first = False
plt.draw()
self.fig.canvas.draw()
self.fig.canvas.flush_events()

View File

@ -27,9 +27,9 @@ import time
import numpy as np
import iqplot
from adq_mr import Adq
from secop.core import Attached, BoolType, Done, FloatRange, HasIodev, \
import secop_psi.iqplot as iqplot
from secop_psi.adq_mr import Adq
from secop.core import Attached, BoolType, Done, FloatRange, HasIO, \
IntRange, Module, Parameter, Readable, StringIO, StringType
from secop.properties import Property
@ -60,7 +60,8 @@ class Roi(Readable):
interval = (0,0)
def initModule(self):
self._main.register_roi(self)
super().initModule()
self.main.register_roi(self)
self.calc_interval()
def calc_interval(self):
@ -90,7 +91,7 @@ class FreqStringIO(StringIO):
end_of_line = '\r'
class Frequency(HasIodev, Readable):
class Frequency(HasIO, Readable):
pars = Attached()
sr = Property('samples per record', datatype=IntRange(), default=16384)
maxy = Property('plot y scale', datatype=FloatRange(), default=0.5)
@ -98,8 +99,9 @@ class Frequency(HasIodev, Readable):
value = Parameter('frequency@I,q', datatype=FloatRange(unit='Hz'), default=0)
basefreq = Parameter('base frequency', FloatRange(unit='Hz'), readonly=False)
nr = Parameter('number of records', datatype=IntRange(1,10000), default=500)
freq = Parameter('target frequency', FloatRange(unit='Hz'), readonly=False, poll=True)
amp = Parameter('amplitude', FloatRange(unit='dBm'), readonly=False, poll=True)
freq = Parameter('target frequency', FloatRange(unit='Hz'), readonly=False)
bw = Parameter('bandwidth lowpassfilter', datatype=FloatRange(unit='Hz'),default=10E6)
amp = Parameter('amplitude', FloatRange(unit='dBm'), readonly=False)
control = Parameter('control loop on?', BoolType(), readonly=False, default=True)
time = Parameter('pulse start time', FloatRange(unit='nsec'),
readonly=False)
@ -116,7 +118,7 @@ class Frequency(HasIodev, Readable):
save = Parameter('save data', BoolType(), readonly=False, default=True)
pollinterval = Parameter(datatype=FloatRange(0,120))
iodevClass = FreqStringIO
ioClass = FreqStringIO
lastfreq = None
old = None
@ -124,9 +126,8 @@ class Frequency(HasIodev, Readable):
interval = (0,0)
def earlyInit(self):
#assert self.iodev.startswith('serial:')
#self._iodev = serial.Serial(self.iodev[7:])
self.adq = Adq(self.nr, self.sr)
super().earlyInit()
self.adq = Adq(self.nr, self.sr, self.bw)
self.roilist = []
self.write_nr(self.nr)
self.skipctrl = 0
@ -155,24 +156,24 @@ class Frequency(HasIodev, Readable):
def set_freq(self):
freq = self.freq + self.basefreq
self.sendRecv('FREQ %.15g;FREQ?' % freq)
self.communicate('FREQ %.15g;FREQ?' % freq)
#self._iodev.readline().decode('ascii')
return freq
def write_amp(self, amp):
reply = self.sendRecv('AMPR %g;AMPR?' % amp)
reply = self.communicate('AMPR %g;AMPR?' % amp)
return float(reply)
def read_amp(self):
reply = self.sendRecv('AMPR?')
reply = self.communicate('AMPR?')
return float(reply)
def write_freq(self, value):
self.skipctrl = 2 # suppress control for the 2 next steps
return value
def read_freq(self):
"""used as main polling loop body"""
def doPoll(self):
"""main poll loop body"""
if self.lastfreq is None:
self.lastfreq = self.set_freq()
self.adq.start()
@ -180,13 +181,19 @@ class Frequency(HasIodev, Readable):
self.starttime = time.time()
times = []
times.append(('init', time.time()))
seadata = {p: float(getattr(self._pars, p)) for p in self._pars.parameters}
seadata = {p: float(getattr(self.pars, p)) for p in self.pars.parameters}
data = self.adq.getdata() # this will wait, if not yet finished
#save sample
#np.save('sample.dat',data)
times.append(('wait',time.time()))
freq = self.lastfreq # data was acquired at this freq
if self.control:
freq = self.lastfreq # data was acquired at this freq
else:
freq = self.set_freq()
seadata['frequency'] = freq
self.lastfreq = self.set_freq()
times.append(('setf',time.time()))
if self.control:
self.lastfreq = self.set_freq()
times.append(('setf',time.time()))
self.adq.start() # start next acq
times.append(('start',time.time()))
roilist = [r for r in self.roilist if r.enable]
@ -194,6 +201,7 @@ class Frequency(HasIodev, Readable):
gates = self.adq.gates_and_curves(data, freq, self.interval,
[r.interval for r in roilist])
if self.save:
times.append(('save',time.time()))
tdata, idata, qdata, pdata = self.adq.curves
seadata['timestep'] = tdata[1] - tdata[0]
iqdata = np.array((idata, qdata, pdata), dtype='f4')
@ -249,4 +257,5 @@ class Frequency(HasIodev, Readable):
self.skipctrl -= 1
elif self.control:
self.freq = sorted((self.freq - self.maxstep, newfreq, self.freq + self.maxstep))[1]
#print(times)
return Done