diff --git a/secop_psi/adq_mr.py b/secop_psi/adq_mr.py new file mode 100644 index 0000000..d0e9fc1 --- /dev/null +++ b/secop_psi/adq_mr.py @@ -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] + diff --git a/secop_psi/dg645.py b/secop_psi/dg645.py index 90e5e24..3c3b413 100644 --- a/secop_psi/dg645.py +++ b/secop_psi/dg645.py @@ -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]) diff --git a/secop_psi/iqplot.py b/secop_psi/iqplot.py new file mode 100644 index 0000000..344f062 --- /dev/null +++ b/secop_psi/iqplot.py @@ -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 .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() diff --git a/secop_psi/ultrasound.py b/secop_psi/ultrasound.py index d239c1f..f645f10 100644 --- a/secop_psi/ultrasound.py +++ b/secop_psi/ultrasound.py @@ -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