WIP new version of ultrasound

Change-Id: Iadb83396a64e277f6f0a37f7a96d92105648c4fe

frappy_psi/adq_mr.py

@@ -1,27 +1,37 @@
-"""
+# *****************************************************************************
-Created on Tue Nov 26 15:42:43 2019
+# This program is free software; you can redistribute it and/or modify it under
-
+# the terms of the GNU General Public License as published by the Free Software
-@author: tartarotti_d-adm
+# Foundation; either version 2 of the License, or (at your option) any later
-"""
+# version.
-
+#
-
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+#
+# You should have received a copy of the GNU General Public License along with
+# this program; if not, write to the Free Software Foundation, Inc.,
+# 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+#
+# Module authors:
+#   Damaris Tartarotti Maimone
+#   Markus Zolliker <markus.zolliker@psi.ch>
+# *****************************************************************************
+"""Wrapper for the ADQ data acquisition card for ultrasound"""
 import sys
 import atexit
 import signal
 import time
 import numpy as np
 import ctypes as ct
-from numpy import sqrt, arctan2, sin, cos
+from scipy.signal import butter, filtfilt
-import scipy.signal
 
-#from pylab import *
-
-#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
+# The following 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_1 = 2
 EXT_TRIG_2 = 7
 EXT_TRIG_3 = 8
 LVL_TRIG = 3
@@ -29,29 +39,36 @@ INT_TRIG = 4
 LVL_FALLING = 0
 LVL_RISING = 1
 
-#samples_per_record=16384
+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
+
 ADQ_TRANSFER_MODE_NORMAL = 0x00
 ADQ_CHANNELS_MASK = 0x3
 
-#f_LO = 40
+GHz = 1e9
-
-def butter_lowpass(cutoff, sr, order=5):
-    nyq = 0.5 * sr
-    normal_cutoff = cutoff / nyq
-    b, a = scipy.signal.butter(order, normal_cutoff, btype = 'low', analog = False)
-    return b, a
 
 
-class Adq(object):
+class Adq:
+    sample_rate = 2 * GHz
     max_number_of_channels = 2
-    samp_freq = 2
+    ndecimate = 50  # decimation ratio (2GHz / 40 MHz)
-    #ndecimate = 50  # decimation ratio (2GHz / 40 MHz)
+    number_of_records = 1
-    ndecimate = 50
+    samples_per_record = 16384
+    bw_cutoff = 10E6
+    trigger = EXT_TRIG_1
+    adq_num = 1
+    UNDEFINED = -1
+    IDLE = 0
+    BUSY = 1
+    READY = 2
+    status = UNDEFINED
+    data = None
+
+    def __init__(self):
+        global ADQAPI
+        ADQAPI = ct.cdll.LoadLibrary("libadq.so.0")
 
-    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
@@ -62,33 +79,29 @@ class Adq(object):
         # 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)
+        n_of_adq = ADQAPI.ADQControlUnit_NofADQ(self.adq_cu)
-        if n_of_ADQ != 1:
+        if n_of_adq != 1:
-            raise ValueError('number of ADQs must be 1, not %d' % n_of_ADQ)
+            raise RuntimeError('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]):
+        if revision[1]:
             print('Local copy')
         else:
             print('SVN Managed')
-            if (revision[2]):
+            if revision[2]:
                 print('Mixed Revision')
             else:
                 print('SVN Updated')
                 print('')
 
-        ADQ_CLOCK_INT_INTREF = 0  #internal clock source
+        ADQAPI.ADQ_SetClockSource(self.adq_cu, self.adq_num, ADQ_CLOCK_EXT_REF)
-        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
@@ -98,38 +111,33 @@ class Adq(object):
         # ADQAPI.ADQ_SetSampleSkip(self.adq_cu, self.adq_num, 1)
         ##########################
 
-        # Set trig mode
+        # set trigger mode
-        self.trigger = EXT_TRIG_1
+        if not ADQAPI.ADQ_SetTriggerMode(self.adq_cu, self.adq_num, self.trigger):
-        #trigger = LVL_TRIG
+            raise RuntimeError('ADQ_SetTriggerMode failed.')
-        success = ADQAPI.ADQ_SetTriggerMode(self.adq_cu, self.adq_num, self.trigger)
+        if self.trigger == LVL_TRIG:
-        if (success == 0):
+            if not ADQAPI.ADQ_SetLvlTrigLevel(self.adq_cu, self.adq_num, -100):
-            print('ADQ_SetTriggerMode failed.')
+                raise RuntimeError('ADQ_SetLvlTrigLevel failed.')
-        if (self.trigger == LVL_TRIG): 
+            if not ADQAPI.ADQ_SetTrigLevelResetValue(self.adq_cu, self.adq_num, 1000):
-            success = ADQAPI.ADQ_SetLvlTrigLevel(self.adq_cu, self.adq_num, -100)
+                raise RuntimeError('ADQ_SetTrigLevelResetValue failed.')
-            if (success == 0):
+            if not ADQAPI.ADQ_SetLvlTrigChannel(self.adq_cu, self.adq_num, 1):
-                print('ADQ_SetLvlTrigLevel failed.')
+                raise RuntimeError('ADQ_SetLvlTrigChannel failed.')
-            success = ADQAPI.ADQ_SetTrigLevelResetValue(self.adq_cu, self.adq_num, 1000)
+            if not ADQAPI.ADQ_SetLvlTrigEdge(self.adq_cu, self.adq_num, LVL_RISING):
-            if (success == 0):
+                raise RuntimeError('ADQ_SetLvlTrigEdge failed.')
-                print('ADQ_SetTrigLevelResetValue failed.')
+        elif self.trigger == EXT_TRIG_1:
-            success = ADQAPI.ADQ_SetLvlTrigChannel(self.adq_cu, self.adq_num, 1)
+            if not ADQAPI.ADQ_SetExternTrigEdge(self.adq_cu, self.adq_num, 2):
-            if (success == 0):
+                raise RuntimeError('ADQ_SetLvlTrigEdge failed.')
-                print('ADQ_SetLvlTrigChannel failed.')
+    #        if not ADQAPI.ADQ_SetTriggerThresholdVoltage(self.adq_cu, self.adq_num, trigger, ct.c_double(0.2)):
-            success = ADQAPI.ADQ_SetLvlTrigEdge(self.adq_cu, self.adq_num, LVL_RISING)
+    #            raise RuntimeError('SetTriggerThresholdVoltage failed.')
-            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()
         atexit.register(self.deletecu)
         signal.signal(signal.SIGTERM, lambda *_: sys.exit(0))
 
-    def setup_target_buffers(self):
+    def init(self, samples_per_record=None, number_of_records=None):
+        """initialize dimensions"""
+        if samples_per_record:
+            self.samples_per_record = samples_per_record
+        if number_of_records:
+            self.number_of_records = number_of_records
         # 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)()
@@ -139,15 +147,11 @@ class Adq(object):
     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);
+        ADQAPI.ADQ_MultiRecordClose(self.adq_cu, self.adq_num)
         # Delete ADQControlunit
         ADQAPI.DeleteADQControlUnit(self.adq_cu)
-        print('ADQ closed')
         
     def start(self):
-        """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,
@@ -155,32 +159,48 @@ class Adq(object):
 
         ADQAPI.ADQ_DisarmTrigger(self.adq_cu, self.adq_num)
         ADQAPI.ADQ_ArmTrigger(self.adq_cu, self.adq_num)
+        self.status = self.BUSY
 
-    def getdata(self):
+    def get_status(self):
-        """wait for aquisition to be finished and get data"""
+        """check if ADQ card is busy"""
-        #start = time.time()
+        if self.status == self.BUSY:
-        while(ADQAPI.ADQ_GetAcquiredAll(self.adq_cu,self.adq_num) == 0):
+            if ADQAPI.ADQ_GetAcquiredAll(self.adq_cu, self.adq_num):
-            time.sleep(0.001)
+                self.status = self.READY
-            #if (self.trigger == SW_TRIG):
+            else:
-            #    ADQAPI.ADQ_SWTrig(self.adq_cu, self.adq_num)
+                if self.trigger == SW_TRIG:
-        #mid = time.time()
+                    ADQAPI.ADQ_SWTrig(self.adq_cu, self.adq_num)
-        status = ADQAPI.ADQ_GetData(self.adq_cu, self.adq_num, self.target_buffers,
+        return self.status
+
+    def get_data(self, dataclass, **kwds):
+        """when ready, get raw data from card, else return cached data
+
+        return
+        """
+        if self.get_status() == self.READY:
+            # Get data from ADQ
+            if not 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);
+                                      0, self.samples_per_record, ADQ_TRANSFER_MODE_NORMAL):
-        #print(time.time()-mid,mid-start)
+                raise RuntimeError('no success from ADQ_GetDATA')
-        if not status:
+            self.data = dataclass(self, **kwds)
-            raise ValueError('no succesS from ADQ_GetDATA')
+            self.status = self.IDLE
-        # Now this is an array with all records, but the time is artificial
+        if self.status == self.UNDEFINED:
+            raise RuntimeError('no data available yet')
+        return self.data
+
+
+class PEdata:
+    def __init__(self, adq):
+        self.sample_rate = adq.sample_rate
+        self.samp_freq = self.sample_rate / GHz
+        self.number_of_records = adq.number_of_records
         data = []
         for ch in range(2):
-            onedim = np.frombuffer(self.target_buffers[ch].contents, dtype=np.int16)
+            onedim = np.frombuffer(adq.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
+            data.append(onedim.reshape(adq.number_of_records, adq.samples_per_record) / float(2**14))  # 14 bits ADC
-        return data
+        # Now this is an array with all records, but the time is artificial
-
+        self.data = data
-    def acquire(self):
-        self.start()
-        return self.getdata()
 
     def sinW(self, sig, freq, ti, tf):
         # sig: signal array
@@ -201,30 +221,28 @@ class Adq(object):
     def mix(self, sigin, sigout, freq, ti, tf):
         # sigin, sigout: signal array, incomping, output
         # freq
-        # ti, tf: initial and end time if sigin
+        # ti, tf: initial and end time of sigin
         a, b = self.sinW(sigin, freq, ti, tf)
-        phase = arctan2(a,b) * 180 / np.pi
+        amp = np.sqrt(a**2 + b**2)
-        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))
+        wave1 = sigout * (a * np.cos(2*np.pi*freq*t) + b * np.sin(2*np.pi*freq*t))
-        wave2 = sigout * (a * sin(2*np.pi*freq*t) - b * cos(2*np.pi*freq*t))
+        wave2 = sigout * (a * np.sin(2*np.pi*freq*t) - b * np.cos(2*np.pi*freq*t))
         return wave1, wave2
 
     def averageiq(self, data, freq, ti, tf):
-        '''Average over records'''
+        """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)
+        # butter lowpass
-
+        nyq = 0.5 * self.sample_rate
-        #ifi = np.array(scipy.signal.filtfilt(b,a,iorq[0]))
+        normal_cutoff = cutoff / nyq
-        #qf = np.array(scipy.signal.filtfilt(b,a,iorq[1]))
+        order = 5
-        iqf = [scipy.signal.filtfilt(b,a,iorq[i]) for i in np.arange(len(iorq))]
+        b, a = butter(order, normal_cutoff, btype='low', analog=False)
-
+        iqf = [filtfilt(b, a, iorq[i]) for i in np.arange(len(iorq))]
         return iqf
 
     def box(self, iorq, ti, tf):
@@ -233,24 +251,63 @@ class Adq(object):
         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):
+    def gates_and_curves(self, freq, pulse, roi, bw_cutoff):
         """return iq values of rois and prepare plottable curves for iq"""
-        self.ndecimate = int(round(2E9/freq))
+        self.ndecimate = int(round(self.sample_rate / freq))
-        times = []
+        # times = []
-        times.append(('aviq', time.time()))
+        # times.append(('aviq', time.time()))
-        iq = self.averageiq(data,freq*1e-9,*pulse)
+        iq = self.averageiq(self.data, freq / GHz, *pulse)
-        times.append(('filtro', time.time()))
+        # times.append(('filtro', time.time()))
-        iqf = self.filtro(iq,self.bw_cutoff)
+        iqf = self.filtro(iq, 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()))
+        # 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])
+        pulsig = np.abs(self.data[0][0])
-        times.append(('pulsig', time.time()))
+        # 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]
 
+
+class RUSdata:
+    def __init__(self, adq, freq, periods):
+        self.sample_rate = adq.sample_rate
+        self.freq = freq
+        self.periods = periods
+        self.samples_per_record = adq.samples_per_record
+        input_signal = np.frombuffer(adq.target_buffers[0].contents, dtype=np.int16)
+        output_signal = np.frombuffer(adq.target_buffers[1].contents, dtype=np.int16)
+        complex_sinusoid = np.exp(1j * 2 * np.pi * self.freq / self.sample_rate * np.arange(len(input_signal)))
+        self.input_mixed = input_signal * complex_sinusoid
+        self.output_mixed = output_signal * complex_sinusoid
+        self.input_mean = self.input_mixed.mean()
+        self.output_mean = self.output_mixed.mean()
+        self.iq = self.output_mean / self.input_mean
+
+    def get_reduced(self, mixed):
+        """get reduced array and normalize"""
+        nper = self.samples_per_record // self.periods
+        mean = mixed.mean()
+        return mixed[:self.period * nper].reshape((-1, nper)).mean(axis=0) / mean
+
+    def calc_quality(self):
+        """get signal quality info
+
+        quality info (small values indicate good quality):
+        - input_std and output_std:
+          the imaginary part indicates deviations in phase
+          the real part indicates deviations in amplitude
+        - input_slope and output_slope:
+          the imaginary part indicates a turning phase (rad/sec)
+          the real part indicates changes in amplitude (0.01 ~= 1%/sec)
+        """
+        reduced = self.get_reduced(self.input_mixed)
+        self.input_stdev = reduced.std()
+
+        reduced = self.get_reduced(self.output_mixed)
+        timeaxis = np.arange(len(reduced)) * self.sample_rate / self.freq
+        self.output_slope = np.polyfit(timeaxis, reduced, 1)[0]

frappy_psi/ultrasound.py

@@ -19,18 +19,19 @@
 """frappy support for ultrasound"""
 
 import math
-#import serial
 import os
 import time
 
 import numpy as np
 
-import frappy_psi.iqplot as iqplot
+from frappy_psi.adq_mr import Adq, PEdata, RUSdata
-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, Writable, Drivable, StringIO, StringType, \
-    IDLE, DISABLED, TupleOf, ArrayOf
+    IDLE, BUSY, DISABLED, ERROR, TupleOf, ArrayOf, Command
 from frappy.properties import Property
+#from frappy.modules import Collector
+
+Collector = Readable
 
 
 def fname_from_time(t, extension):
@@ -92,80 +93,22 @@ class FreqStringIO(StringIO):
     end_of_line = '\r'
 
 
-class Frequency(HasIO, Readable):
+class Frequency(HasIO, Writable):
-    pars = Attached()
+    value = Parameter('frequency', unit='Hz')
-    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'),
-                    readonly=False)
-    pulselen = Parameter('adjusted pulse length (integer number of periods)', FloatRange(unit='nsec'), default=1)
-    maxstep = Parameter('max frequency step', FloatRange(unit='Hz'), readonly=False,
-                      default=10000)
-    minstep = Parameter('min frequency step for slope calculation', FloatRange(unit='Hz'),
-                        readonly=False, default=4000)
-    slope = Parameter('inphase/frequency slope', FloatRange(), readonly=False,
-                      default=1e6)
-    plot = Parameter('create plot images', BoolType(), readonly=False, default=True)
-    save = Parameter('save data', BoolType(), readonly=False, default=True)
-    pollinterval = Parameter(datatype=FloatRange(0,120))
 
+    last_change = 0
     ioClass = FreqStringIO
+    dif = None
 
-    lastfreq = None
+    def register_dif(self, dif):
-    old = None
+        self.dif = dif
-    starttime = None
-    interval = (0,0)
 
-    def earlyInit(self):
+    def write_target(self, value):
-        super().earlyInit()
+        self.communicate('FREQ %.15g;FREQ?' % value)
-        self.adq = Adq(self.nr, self.sr, self.bw)
+        self.last_change = time.time()
-        self.roilist = []
+        if self.dif:
-        self.write_nr(self.nr)
+            self.dif.read_value()
-        self.write_sr(self.sr)
-        self.skipctrl = 0
-        self.plotter = iqplot.Plot(self.maxy)
-        self.calc_interval()
-
-    def calc_interval(self):
-        self.interval = (self.time, self.time + self.size)
-
-    def write_time(self, value):
-        self.time = value
-        self.calc_interval()
-        return Done
-
-    def write_size(self, value):
-        self.size = value
-        self.calc_interval()
-        return Done
-
-    def write_nr(self, value):
-        # self.pollinterval = value * 0.0001
-        return value
-
-    def write_sr(self, value):
-        return value
-
-    def register_roi(self, roi):
-        self.roilist.append(roi)
-
-    def set_freq(self):
-        freq = self.freq + self.basefreq
-        self.communicate('FREQ %.15g;FREQ?' % freq)
-        #self._iodev.readline().decode('ascii')
-        return freq
 
     def write_amp(self, amp):
         reply = self.communicate('AMPR %g;AMPR?' % amp)
@@ -175,105 +118,196 @@ class Frequency(HasIO, Readable):
         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 doPoll(self):
+class FrequencyDif(Readable):
-        """main poll loop body"""
+    freq = Attached(Frequency)
-        if self.lastfreq is None:
+    base = Parameter('base frequency', FloatRange(unit='Hz'), default=0)
-            self.lastfreq = self.set_freq()
+    value = Parameter('difference to base frequency', FloatRange(unit='Hz'), default=0)
-            if self.rusmode:
+
-                self.sr = int(12e9/self.lastfreq) #picking up 12 period at the ith frequency in the time scale
+    def initModule(self):
-#                self.adq.samples_per_record = self.sr
+        super().initModule()
+        self.freq.register_dif(self)
+
+    def read_value(self):
+        return self.freq - self.base
+
+
+class Base(Collector):
+    freq = Attached()
+    adq = Attached(Adq)
+    sr = Parameter('samples per record', datatype=IntRange(1, 1E9), default=16384)
+    pollinterval = Parameter(datatype=FloatRange(0, 120))  # allow pollinterval = 0
+    _data = None
+    _data_args = None
+
+    def read_status(self):
+        adqstate = self.adq.get_status()
+        if adqstate == Adq.BUSY:
+            return BUSY, 'acquiring'
+        if adqstate == Adq.UNDEFINED:
+            return ERROR, 'no data yet'
+        if adqstate == Adq.READY:
+            return IDLE, 'new data available'
+        return IDLE, ''
+
+    def get_data(self):
+        data = self.adq.get_data(*self._data_args)
+        if id(data) != id(self._data):
+            self._data = data
+            return data
+        return None
+
+
+class PulseEcho(Base):
+    value = Parameter("t, i, q, pulse curves",
+                      TupleOf(*[ArrayOf(FloatRange(), 0, 16283) for _ in range(4)]), default=[[]] * 4)
+    nr = Parameter('number of records', datatype=IntRange(1, 9999), default=500)
+    bw = Parameter('bandwidth lowpassfilter', datatype=FloatRange(unit='Hz'), default=10E6)
+    control = Parameter('control loop on?', BoolType(), readonly=False, default=True)
+    time = Parameter('pulse start time', FloatRange(unit='nsec'),
+                     readonly=False)
+    size = Parameter('pulse length (starting from time)', FloatRange(unit='nsec'),
+                     readonly=False)
+    pulselen = Parameter('adjusted pulse length (integer number of periods)', FloatRange(unit='nsec'), default=1)
+
+    starttime = None
+
+    def initModule(self):
+        super().initModule()
+        self.adq = Adq()
+        self.adq.init(self.sr, self.nr)
+        self.roilist = []
+
+    def write_nr(self, value):
+        self.adq.init(self.sr, value)
+
+    def write_sr(self, value):
+        self.adq.init(value, self.nr)
+
+    def write_bw(self, value):
+        self.adq.bw_cutoff = value
+
+    def register_roi(self, roi):
+        self.roilist.append(roi)
+
+    def go(self):
+        self.starttime = time.time()
         self.adq.start()
 
-        if self.starttime is None:
+    def read_value(self):
-            self.starttime = time.time()
+        if self.get_rawdata():  # new data available
-        times = []
-        times.append(('init', time.time()))
-        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()))
-        if self.control:
-            freq = self.lastfreq  # data was acquired at this freq
-        else:
-            freq = self.set_freq()
-        seadata['frequency'] = freq
-        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]
-
+            freq = self.freq.value
-        gates = self.adq.gates_and_curves(data, freq, self.interval,
+            gates = self.adq.gates_and_curves(self._data, freq,
+                                              (self.time, self.time + self.size),
                                               [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
-            seadata['timestep'] = tdata[1] - tdata[0]
-            iqdata = np.array((idata, qdata, pdata), dtype='f4')
-            ts = seadata['timestamp']
-            if ts:
-                filename = fname_from_time(ts, '.npz')
-                seanp = np.array(list(seadata.items()), dtype=[('name', 'U16'), ('value', 'f8')])
-                np.savez(filename, seadata=seanp, iqdata=iqdata)
-                # can be load back via
-                # content = np.load(filename)
-                # content['seadata'], content['iqdata']
-        self.pulselen = self.adq.pulselen
-        times.append(('ana',time.time()))
-        if self.plot:
-            # get reduced interval from adq.sinW
-            pulseint = (self.interval[0], self.interval[0] + self.pulselen)
-            try:
-                self.plotter.plot(
-                    self.adq.curves,
-                    rois=[pulseint] + [r.interval for r in roilist],
-                    average=([r.time for r in roilist],
-                             [r.i for r in roilist],
-                             [r.q for r in roilist]))
-            except Exception as e:
-                self.log.warning('can not plot %r' % e)
-        else:
-            self.plotter.close()
-        now = time.time()
-        times.append(('plot',now))
-        # print(' '.join('%s %5.3f' % (label, t - self.starttime) for label, t in times))
-        self.starttime = now
-        self.value = freq - self.basefreq
             for i, roi in enumerate(roilist):
                 roi.i = a = gates[i][0]
                 roi.q = b = gates[i][1]
                 roi.value = math.sqrt(a ** 2 + b ** 2)
                 roi.phase = math.atan2(a, b) * 180 / math.pi
-        inphase = self.roilist[0].i
+        return self.adq.curves
-        if self.control:
+
-            newfreq = freq + inphase * self.slope - self.basefreq
+        # TODO: CONTROL
-            # step = sorted((-self.maxstep, inphase * self.slope, self.maxstep))[1]
+        # inphase = self.roilist[0].i
-        if self.old:
+        # if self.control:
-            fdif = freq - self.old[0]
+        #     newfreq = freq + inphase * self.slope - self.basefreq
-            idif = inphase - self.old[1]
+        #     # step = sorted((-self.maxstep, inphase * self.slope, self.maxstep))[1]
-            if abs(fdif) >= self.minstep:
+        # if self.old:
-                self.slope = - fdif / idif
+        #     fdif = freq - self.old[0]
-        else:
+        #     idif = inphase - self.old[1]
-            fdif = 0
+        #     if abs(fdif) >= self.minstep:
-            idif = 0
+        #         self.slope = - fdif / idif
-            newfreq = freq + self.minstep
+        # else:
-        self.old = (freq, inphase)
+        #     fdif = 0
-        if self.skipctrl > 0:  # do no control for some time after changing frequency
+        #     idif = 0
-            self.skipctrl -= 1
+        #     newfreq = freq + self.minstep
-        elif self.control:
+        # self.old = (freq, inphase)
-            self.freq = sorted((self.freq - self.maxstep, newfreq, self.freq + self.maxstep))[1]
+        # if self.skipctrl > 0:  # do no control for some time after changing frequency
-        #print(times)
+        #     self.skipctrl -= 1
-        return Done
+        # elif self.control:
+        #     self.freq = sorted((self.freq - self.maxstep, newfreq, self.freq + self.maxstep))[1]
 
 
-class Curves(Readable):
+class RUS(Base):
-    value = Parameter("t, i, q, pulse curves",
+    value = Parameter('averaged (I, Q) tuple', TupleOf(FloatRange(), FloatRange()))
-                      TupleOf(*[ArrayOf(FloatRange(), 0, 16283) for _ in range(4)]), default=[[]] * 4)
+    periods = Parameter('number of periods', IntRange(1, 9999), default=12)
+    scale = Parameter('scale,taking into account input attenuation', FloatRange(), default=0.1)
+    input_phase_stddev = Parameter('input signal quality', FloatRange(unit='rad'))
+    output_phase_slope = Parameter('output signal phase slope', FloatRange(unit='rad/sec'))
+    output_amp_slope = Parameter('output signal amplitude change', FloatRange(unit='1/sec'))
+    phase = Parameter('phase', FloatRange(unit='deg'))
+    amp = Parameter('amplitude', FloatRange())
+
+    starttime = None
+    _data_args = None
+
+    def initModule(self):
+        super().initModule()
+        self.adq = Adq()
+        # self.write_periods(self.periods)
+
+    def read_value(self):
+        if self._data_args is None:
+            return self.value  # or may we raise as no value is defined yet?
+        data = self.get_data(RUSdata, *self._data_args)
+        if data:
+            # data available
+            data.calc_quality()
+            self.input_phase_stddev = data.input_stddev.imag
+            self.output_phase_slope = data.output_slope.imag
+            self.output_amp_slope = data.output_slope.real
+
+            iq = data.iq * self.scale
+            self.phase = np.arctan2(iq.imag, iq.real) * 180 / np.pi
+            self.amp = np.abs(iq.imag, iq.real)
+            return iq.real, iq.imag
+        return self.value
+
+    def go(self):
+        self.starttime = time.time()
+        freq = self.freq.value
+        self._data_args = (RUSdata, freq, self.periods)
+        self.sr = round(self.periods * self.adq.sample_rate / freq)
+        self.adq.init(self.sr, 1)
+        self.adq.start()
+        self.read_status()
+
+
+class ControlLoop:
+    maxstep = Parameter('max frequency step', FloatRange(unit='Hz'), readonly=False,
+                      default=10000)
+    minstep = Parameter('min frequency step for slope calculation', FloatRange(unit='Hz'),
+                        readonly=False, default=4000)
+    slope = Parameter('inphase/frequency slope', FloatRange(), readonly=False,
+                      default=1e6)
+
+
+# class Frequency(HasIO, Readable):
+#     pars = Attached()
+#     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'),
+#                     readonly=False)
+#     pulselen = Parameter('adjusted pulse length (integer number of periods)', FloatRange(unit='nsec'), default=1)
+#     maxstep = Parameter('max frequency step', FloatRange(unit='Hz'), readonly=False,
+#                       default=10000)
+#     minstep = Parameter('min frequency step for slope calculation', FloatRange(unit='Hz'),
+#                         readonly=False, default=4000)
+#     slope = Parameter('inphase/frequency slope', FloatRange(), readonly=False,
+#                       default=1e6)
+#     plot = Parameter('create plot images', BoolType(), readonly=False, default=True)
+#     save = Parameter('save data', BoolType(), readonly=False, default=True)
+#     pollinterval = Parameter(datatype=FloatRange(0,120))