diff --git a/cfg/PEUS.py b/cfg/PEUS.py new file mode 100644 index 0000000..164276d --- /dev/null +++ b/cfg/PEUS.py @@ -0,0 +1,67 @@ +Node(equipment_id = 'pe_ultrasound.psi.ch', + description = 'pulse echo ultra sound setup', + interface = 'tcp://5000', + ) + +Mod('f', + cls = 'frappy_psi.ultrasound.Frequency', + description = 'ultrasound frequency and acquisition loop', + uri = 'serial:///dev/ttyS1', + pars = 'pars', + pollinterval = 0.1, + time = 900, # start time + size = 5000, + freq = 1.17568e+06, + basefreq = 4.14902e+07, + control = False, + rusmode = False, + amp = 5.0, + nr = 1000, #500 #300 #100 #50 #30 #10 #5 #3 #1 #1000 #500 #300 #100 #50 #30 #10 #5 #3 #1 #500 + sr = 32768, #16384 + plot = True, + maxstep = 100000, + bw = 10E6, #butter worth filter bandwidth + maxy = 0.7, # y scale for plot + curves = 'curves', # module to transmit curves: + ) + +Mod('curves', + cls = 'frappy_psi.ultrasound.Curves', + description = 't, I, Q and pulse arrays for plot', + ) + +Mod('delay', + cls = 'frappy__psi.dg645.Delay', + description = 'delay line with 2 channels', + uri = 'serial:///dev/ttyS2', + on1 = 1e-9, + on2 = 1E-9, + off1 = 400e-9, + off2 = 600e-9, + ) + +Mod('pars', + cls = 'frappy_psi.ultrasound.Pars', + description = 'SEA parameters', + ) + +def roi(nr, time=None, size=300): + Mod(f'roi{nr}', + cls = 'frappy_psi.ultrasound.Roi', + description = f'I/Q of region {nr}', + main = 'f', + time=time or 4000, + size=size, + enable=time is not None, + ) + +roi(0, 2450) # you may add size as argument if not default +roi(1, 5950) +roi(2, 9475) +roi(3, 12900) +roi(4, 16100) +roi(5) # disabled +roi(6) +roi(7) +roi(8) +roi(9) diff --git a/cfg/RUS.py b/cfg/RUS.py new file mode 100644 index 0000000..9b5f235 --- /dev/null +++ b/cfg/RUS.py @@ -0,0 +1,62 @@ +Node(equipment_id = 'r_ultrasound.psi.ch', + description = 'resonant ultra sound setup', + interface = 'tcp://5000', +) + +Mod('f', + cls = 'frappy_psi.ultrasound.Frequency', + description = 'ultrasound frequency and acquisition loop', + uri = 'serial:///dev/ttyS1', + pars = 'pars', + pollinterval = 0.1, + time = 900, # start time + size = 5000, + freq = 1.e+03, + basefreq = 1.E+3, + control = False, + rusmode = False, + amp = 2.5, + nr = 1, #500 #300 #100 #50 #30 #10 #5 #3 #1 #1000 #500 #300 #100 #50 #30 #10 #5 #3 #1 #500 + sr = 1E8, #16384 + plot = True, + maxstep = 100000, + bw = 10E6, #butter worth filter bandwidth + maxy = 0.7, # y scale for plot + curves = 'curves', # module to transmit curves: + ) + +Mod('curves', + cls = 'frappy_psi.ultrasound.Curves', + description = 't, I, Q and pulse arrays for plot', + ) + +Mod('roi0', + cls = 'frappy_psi.ultrasound.Roi', + description = 'I/Q of region in the control loop', + time = 300, # this is the center of roi: + size = 5000, + main = f, + ) + +Mod('roi1', + cls = 'frappy_psi.ultrasound.Roi', + description = 'I/Q of region 1', + time = 100, # this is the center of roi: + size = 300, + main = f, + ) + +Mod('delay', + cls = 'frappy__psi.dg645.Delay', + description = 'delay line with 2 channels', + uri = 'serial:///dev/ttyS2', + on1 = 1e-9, + on2 = 1E-9, + off1 = 400e-9, + off2 = 600e-9, + ) + +Mod('pars', + cls = 'frappy_psi.ultrasound.Pars', + description = 'SEA parameters', + ) diff --git a/frappy_psi/adq_mr.py b/frappy_psi/adq_mr.py index 0786ba4..1f85b11 100644 --- a/frappy_psi/adq_mr.py +++ b/frappy_psi/adq_mr.py @@ -5,15 +5,17 @@ Created on Tue Nov 26 15:42:43 2019 """ +import sys +import atexit +import signal +import time import numpy as np import ctypes as ct -import time from numpy import sqrt, arctan2, sin, cos +import scipy.signal #from pylab import * -from scipy import signal - #ADQAPI = ct.cdll.LoadLibrary("ADQAPI.dll") ADQAPI = ct.cdll.LoadLibrary("libadq.so.0") @@ -36,7 +38,7 @@ ADQ_CHANNELS_MASK = 0x3 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) + b, a = scipy.signal.butter(order, normal_cutoff, btype = 'low', analog = False) return b, a @@ -82,12 +84,12 @@ class Adq(object): 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) @@ -95,7 +97,7 @@ class Adq(object): # Sample skip #ADQAPI.ADQ_SetSampleSkip(self.adq_cu, self.adq_num, 1) ########################## - + # Set trig mode self.trigger = EXT_TRIG_1 #trigger = LVL_TRIG @@ -105,13 +107,13 @@ class Adq(object): if (self.trigger == LVL_TRIG): success = ADQAPI.ADQ_SetLvlTrigLevel(self.adq_cu, self.adq_num, -100) if (success == 0): - print('ADQ_SetLvlTrigLevel failed.') + print('ADQ_SetLvlTrigLevel failed.') success = ADQAPI.ADQ_SetTrigLevelResetValue(self.adq_cu, self.adq_num, 1000) if (success == 0): - print('ADQ_SetTrigLevelResetValue failed.') + print('ADQ_SetTrigLevelResetValue failed.') success = ADQAPI.ADQ_SetLvlTrigChannel(self.adq_cu, self.adq_num, 1) if (success == 0): - print('ADQ_SetLvlTrigChannel failed.') + print('ADQ_SetLvlTrigChannel failed.') success = ADQAPI.ADQ_SetLvlTrigEdge(self.adq_cu, self.adq_num, LVL_RISING) if (success == 0): print('ADQ_SetLvlTrigEdge failed.') @@ -124,35 +126,38 @@ class Adq(object): # print('SetTriggerThresholdVoltage failed.') print("CHANNEL:"+str(ct.c_int(ADQAPI.ADQ_GetLvlTrigChannel(self.adq_cu, self.adq_num)))) self.setup_target_buffers() - + atexit.register(self.deletecu) + signal.signal(signal.SIGTERM, lambda *_: sys.exit(0)) + 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) - + print('ADQ closed') + def start(self): - """start datat acquisition""" + """start data 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""" + """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) @@ -163,7 +168,7 @@ class Adq(object): 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) + #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 @@ -176,52 +181,7 @@ class Adq(object): 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 @@ -257,16 +217,16 @@ class Adq(object): 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))] + #ifi = np.array(scipy.signal.filtfilt(b,a,iorq[0])) + #qf = np.array(scipy.signal.filtfilt(b,a,iorq[1])) + iqf = [scipy.signal.filtfilt(b,a,iorq[i]) for i in np.arange(len(iorq))] + + return iqf - return iqf - def box(self, iorq, ti, tf): si = int(self.samp_freq * ti) sf = int(self.samp_freq * tf) @@ -275,12 +235,15 @@ class Adq(object): def gates_and_curves(self, data, freq, pulse, roi): """return iq values of rois and prepare plottable curves for iq""" + self.ndecimate = int(round(2E9/freq)) 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 + ll = m * self.ndecimate + iqf = [iqfx[0:ll] for iqfx in iqf] 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 diff --git a/frappy_psi/ultrasound.py b/frappy_psi/ultrasound.py index b2e0cdf..a35559e 100644 --- a/frappy_psi/ultrasound.py +++ b/frappy_psi/ultrasound.py @@ -28,7 +28,8 @@ import numpy as np import frappy_psi.iqplot as iqplot from frappy_psi.adq_mr import Adq from frappy.core import Attached, BoolType, Done, FloatRange, HasIO, \ - IntRange, Module, Parameter, Readable, StringIO, StringType + IntRange, Module, Parameter, Readable, StringIO, StringType, \ + IDLE, DISABLED, TupleOf, ArrayOf from frappy.properties import Property @@ -52,7 +53,6 @@ class Roi(Readable): time = Parameter('start time', FloatRange(unit='nsec'), readonly=False) size = Parameter('interval (symmetric around time)', FloatRange(unit='nsec'), readonly=False) enable = Parameter('calculate this roi', BoolType(), readonly=False, default=True) - #status = Parameter(export=False) pollinterval = Parameter(export=False) interval = (0,0) @@ -65,6 +65,9 @@ class Roi(Readable): def calc_interval(self): self.interval = (self.time - 0.5 * self.size, self.time + 0.5 * self.size) + def read_status(self): + return (IDLE, '') if self.enable else (DISABLED, 'disabled') + def write_time(self, value): self.time = value self.calc_interval() @@ -82,7 +85,7 @@ class Pars(Module): timestamp = Parameter('unix timestamp', StringType(), default='0', readonly=False) temperature = Parameter('T', FloatRange(unit='K'), default=0, readonly=False) mf = Parameter('field', FloatRange(unit='T'), default=0, readonly=False) - sr = Parameter('rotaion angle', FloatRange(unit='deg'), default=0, readonly=False) + sr = Parameter('rotation angle', FloatRange(unit='deg'), default=0, readonly=False) class FreqStringIO(StringIO): @@ -91,16 +94,18 @@ class FreqStringIO(StringIO): class Frequency(HasIO, Readable): pars = Attached() - sr = Property('samples per record', datatype=IntRange(), default=16384) + curves = Attached(mandatory=False) maxy = Property('plot y scale', datatype=FloatRange(), default=0.5) 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) + sr = Parameter('samples per record', datatype=IntRange(1,1E9), default=16384) 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) + rusmode = Parameter('RUS mode on?', BoolType(), readonly=False, default=False) time = Parameter('pulse start time', FloatRange(unit='nsec'), readonly=False) size = Parameter('pulse length (starting from time)', FloatRange(unit='nsec'), @@ -128,6 +133,7 @@ class Frequency(HasIO, Readable): self.adq = Adq(self.nr, self.sr, self.bw) self.roilist = [] self.write_nr(self.nr) + self.write_sr(self.sr) self.skipctrl = 0 self.plotter = iqplot.Plot(self.maxy) self.calc_interval() @@ -149,6 +155,9 @@ class Frequency(HasIO, Readable): # self.pollinterval = value * 0.0001 return value + def write_sr(self, value): + return value + def register_roi(self, roi): self.roilist.append(roi) @@ -174,7 +183,11 @@ class Frequency(HasIO, Readable): """main poll loop body""" if self.lastfreq is None: self.lastfreq = self.set_freq() + if self.rusmode: + self.sr = int(12e9/self.lastfreq) #picking up 12 period at the ith frequency in the time scale +# self.adq.samples_per_record = self.sr self.adq.start() + if self.starttime is None: self.starttime = time.time() times = [] @@ -198,6 +211,8 @@ class Frequency(HasIO, Readable): gates = self.adq.gates_and_curves(data, freq, self.interval, [r.interval for r in roilist]) + if self.curves: # if attached Curves module is defined, update it + self.curves.value = self.adq.curves if self.save: times.append(('save',time.time())) tdata, idata, qdata, pdata = self.adq.curves @@ -257,3 +272,8 @@ class Frequency(HasIO, Readable): self.freq = sorted((self.freq - self.maxstep, newfreq, self.freq + self.maxstep))[1] #print(times) return Done + + +class Curves(Readable): + value = Parameter("t, i, q, pulse curves", + TupleOf(*[ArrayOf(FloatRange(), 0, 16283) for _ in range(4)]), default=[[]] * 4)