From f8c52af3ac1a59b78aec9ddc8ec1b704c4ced43d Mon Sep 17 00:00:00 2001 From: Markus Zolliker Date: Mon, 17 Mar 2025 09:37:13 +0100 Subject: [PATCH] frappy_psi.ultrasound: after rework (still wip) Change-Id: I200cbeca2dd0f030a01a78ba4d38c342c3c8c8e3 --- frappy_psi/adq_mr.py | 144 ++++++++++++++++------------ frappy_psi/ultrasound.py | 202 +++++++++++++++++++++++++++------------ 2 files changed, 224 insertions(+), 122 deletions(-) diff --git a/frappy_psi/adq_mr.py b/frappy_psi/adq_mr.py index 787ba9c..f40e858 100644 --- a/frappy_psi/adq_mr.py +++ b/frappy_psi/adq_mr.py @@ -58,12 +58,8 @@ class Adq: bw_cutoff = 10E6 trigger = EXT_TRIG_1 adq_num = 1 - UNDEFINED = -1 - IDLE = 0 - BUSY = 1 - READY = 2 - status = UNDEFINED data = None + busy = False def __init__(self): global ADQAPI @@ -84,7 +80,11 @@ class Adq: ADQAPI.ADQControlUnit_FindDevices(self.adq_cu) n_of_adq = ADQAPI.ADQControlUnit_NofADQ(self.adq_cu) if n_of_adq != 1: - raise RuntimeError('number of ADQs must be 1, not %d' % n_of_adq) + print('number of ADQs must be 1, not %d' % n_of_adq) + print('it seems the ADQ was not properly closed') + raise RuntimeError('please try again or reboot') + atexit.register(self.deletecu) + signal.signal(signal.SIGTERM, lambda *_: sys.exit(0)) rev = ADQAPI.ADQ_GetRevision(self.adq_cu, self.adq_num) revision = ct.cast(rev, ct.POINTER(ct.c_int)) @@ -129,11 +129,9 @@ class Adq: # if not ADQAPI.ADQ_SetTriggerThresholdVoltage(self.adq_cu, self.adq_num, trigger, ct.c_double(0.2)): # raise RuntimeError('SetTriggerThresholdVoltage failed.') print("CHANNEL:"+str(ct.c_int(ADQAPI.ADQ_GetLvlTrigChannel(self.adq_cu, self.adq_num)))) - atexit.register(self.deletecu) - signal.signal(signal.SIGTERM, lambda *_: sys.exit(0)) def init(self, samples_per_record=None, number_of_records=None): - """initialize dimensions""" + """initialize dimensions and store result object""" if samples_per_record: self.samples_per_record = samples_per_record if number_of_records: @@ -145,13 +143,20 @@ class Adq: bufp.contents = (ct.c_int16 * self.samples_per_record * self.number_of_records)() def deletecu(self): + cu = self.__dict__.pop('adq_cu', None) + if cu is None: + print('deletecu already called') + return + print('deletecu called') # 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_DisarmTrigger(cu, self.adq_num) + ADQAPI.ADQ_MultiRecordClose(cu, self.adq_num) + print('delete cu') # Delete ADQControlunit - ADQAPI.DeleteADQControlUnit(self.adq_cu) - - def start(self): + ADQAPI.DeleteADQControlUnit(cu) + print('cu deleted') + + def start(self, data): # Start acquisition ADQAPI.ADQ_MultiRecordSetup(self.adq_cu, self.adq_num, self.number_of_records, @@ -159,35 +164,34 @@ class Adq: ADQAPI.ADQ_DisarmTrigger(self.adq_cu, self.adq_num) ADQAPI.ADQ_ArmTrigger(self.adq_cu, self.adq_num) - self.status = self.BUSY + self.data = data - def get_status(self): - """check if ADQ card is busy""" - if self.status == self.BUSY: - if ADQAPI.ADQ_GetAcquiredAll(self.adq_cu, self.adq_num): - self.status = self.READY - else: - if self.trigger == SW_TRIG: - ADQAPI.ADQ_SWTrig(self.adq_cu, self.adq_num) - return self.status + def get_data(self): + """get new data if available""" + ready = False + data = self.data + if not data: + self.busy = False + return None # no new data - 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): - raise RuntimeError('no success from ADQ_GetDATA') - self.data = dataclass(self, **kwds) - self.status = self.IDLE - if self.status == self.UNDEFINED: - raise RuntimeError('no data available yet') - return self.data + if ADQAPI.ADQ_GetAcquiredAll(self.adq_cu, self.adq_num): + ready = True + else: + if self.trigger == SW_TRIG: + ADQAPI.ADQ_SWTrig(self.adq_cu, self.adq_num) + if not ready: + self.busy = True + return None + # 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): + raise RuntimeError('no success from ADQ_GetDATA') + self.data = None + data.retrieve(self) + return data class PEdata: @@ -230,7 +234,6 @@ class PEdata: wave1 = sigout * (a * np.cos(2*np.pi*freq*t) + b * np.sin(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""" iorq = np.array([self.mix(data[0][i], data[1][i], freq, ti, tf) for i in range(self.number_of_records)]) @@ -273,41 +276,56 @@ class PEdata: return [self.box(iqf, *r) for r in roi] +class Namespace: + """holds channel or other data""" + def __init__(self, **kwds): + self.__dict__.update(**kwds) + + 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 + self.inp = Namespace(idx=0, name='input') + self.out = Namespace(idx=1, name='output') + self.channels = (self.inp, self.out) - 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 retrieve(self, adq): + npts = self.samples_per_record + complex_sinusoid = np.exp(1j * 2 * np.pi * self.freq / self.sample_rate * np.arange(npts)) + for chan in self.channels: # looping over input and output + signal = np.frombuffer(adq.target_buffers[chan.idx].contents, dtype=np.int16) + chan.ovr_rate = (np.abs(signal) > 8000).sum() / npts # fraction of points outside range + # calculate RMS * sqrt(2) -> peak sinus amplitude. + # may be higher than the input range by a factor 1.4 when heavily clipped + # we need to convert signal from int16 to float -> use 2.0 as exponent does the trick + chan.amplitude = np.sqrt((signal ** 2.0).mean()) + chan.mixed = signal * complex_sinusoid + chan.mean = chan.mixed.mean() + if self.inp.mean: + self.iq = self.out.mean / self.inp.mean + else: + self.iq = 0 - def calc_quality(self): + def get_quality(self): """get signal quality info quality info (small values indicate good quality): - - input_std and output_std: + - input_stddev: the imaginary part indicates deviations in phase the real part indicates deviations in amplitude - - input_slope and output_slope: + - 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() + nper = self.samples_per_record // self.periods + for chan in self.channels: + mean = chan.mixed.mean() + chan.reduced = chan.mixed[:self.periods * nper].reshape((-1, nper)).mean(axis=0) / mean - 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] + timeaxis = np.arange(len(self.out.reduced)) * self.sample_rate / self.freq + return Namespace( + input_stddev = self.inp.reduced.std(), + output_slope = np.polyfit(timeaxis, self.out.reduced, 1)[0]) diff --git a/frappy_psi/ultrasound.py b/frappy_psi/ultrasound.py index ec46528..e817c2c 100644 --- a/frappy_psi/ultrasound.py +++ b/frappy_psi/ultrasound.py @@ -25,9 +25,9 @@ import time import numpy as np from frappy_psi.adq_mr import Adq, PEdata, RUSdata -from frappy.core import Attached, BoolType, Done, FloatRange, HasIO, \ +from frappy.core import Attached, BoolType, Done, FloatRange, HasIO, StatusType, \ IntRange, Module, Parameter, Readable, Writable, Drivable, StringIO, StringType, \ - IDLE, BUSY, DISABLED, ERROR, TupleOf, ArrayOf, Command + IDLE, BUSY, DISABLED, WARN, ERROR, TupleOf, ArrayOf, Command, Attached from frappy.properties import Property #from frappy.modules import Collector @@ -100,15 +100,22 @@ class Frequency(HasIO, Writable): last_change = 0 ioClass = FreqStringIO dif = None + _freq = None def register_dif(self, dif): self.dif = dif + def read_value(self): + if self._freq is None: + self._freq = float(self.communicate('FREQ?')) + return self._freq + def write_target(self, value): - self.communicate('FREQ %.15g;FREQ?' % value) + self._freq = float(self.communicate('FREQ %.15g;FREQ?' % value)) self.last_change = time.time() if self.dif: self.dif.read_value() + return self._freq def write_amp(self, amp): reply = self.communicate('AMPR %g;AMPR?' % amp) @@ -132,30 +139,18 @@ class FrequencyDif(Readable): return self.freq - self.base -class Base(Collector): +class Base: 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 + adq = 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 shutdownModule(self): + if self.adq: + print('shutdownModule') + self.adq.deletecu() + print('shutdoneModule done') + self.adq = None - 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): @@ -170,7 +165,7 @@ class PulseEcho(Base): readonly=False) pulselen = Parameter('adjusted pulse length (integer number of periods)', FloatRange(unit='nsec'), default=1) - starttime = None + _starttime = None def initModule(self): super().initModule() @@ -190,11 +185,13 @@ class PulseEcho(Base): def register_roi(self, roi): self.roilist.append(roi) - def go(self): - self.starttime = time.time() - self.adq.start() + # TODO: fix + # def go(self): + # self._starttime = time.time() + # self.adq.start() def read_value(self): + # TODO: data = self.get_data() if self.get_rawdata(): # new data available roilist = [r for r in self.roilist if r.enable] freq = self.freq.value @@ -229,54 +226,141 @@ class PulseEcho(Base): # self.freq = sorted((self.freq - self.maxstep, newfreq, self.freq + self.maxstep))[1] -class RUS(Base): - value = Parameter('averaged (I, Q) tuple', TupleOf(FloatRange(), FloatRange())) - 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()) +CONTINUE = 0 +GO = 1 +DONE_GO = 2 +WAIT_GO = 3 - starttime = None - _data_args = None + +class RUS(Base, Collector): + freq = Attached() + imod = Attached(mandatory=False) + qmod = Attached(mandatory=False) + value = Parameter('averaged (I, Q) tuple', TupleOf(FloatRange(), FloatRange())) + status = Parameter(datatype=StatusType(Readable, 'BUSY')) + periods = Parameter('number of periods', IntRange(1, 9999), default=12) + input_range = Parameter('input range (taking in to account attenuation)', FloatRange(unit='V'), default=0.1, readonly=False) + output_range = Parameter('output range', FloatRange(unit='V'), default=0.1, readonly=False) + input_phase_stddev = Parameter('input signal quality', FloatRange(unit='rad'), default=0) + output_phase_slope = Parameter('output signal phase slope', FloatRange(unit='rad/sec'), default=0) + output_amp_slope = Parameter('output signal amplitude change', FloatRange(unit='1/sec'), default=0) + input_amplitude = Parameter('input signal amplitude', FloatRange(unit='V'), default=0) + output_amplitude = Parameter('output signal amplitude', FloatRange(unit='V'), default=0) + phase = Parameter('phase', FloatRange(unit='deg'), default=0) + amp = Parameter('amplitude', FloatRange(), default=0) + continuous = Parameter('continuous mode', BoolType(), readonly=False, default=True) + pollinterval = Parameter(default=1) + + _starttime = None + _iq = 0 + _wait_until = 0 # deadline for returning to continuous mode + _action = CONTINUE # one of CONTINUE, GO, DONE_GO, WAIT_GO + _status = IDLE, 'no data yet' + _busy = False # waiting for end of aquisition (not the same as self.status[0] == BUSY) def initModule(self): super().initModule() self.adq = Adq() + self._ovr_rate = {} # 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 + def doPoll(self): + try: + data = self.adq.get_data() + except Exception as e: + self.set_status(ERROR, repr(e)) + self._busy = False + self._action = WAIT_GO + self.wait_until = time.time() + 2 + return - iq = data.iq * self.scale + self.setFastPoll(False) + if data: # this is new data + self._data = data + for chan in data.channels: + if chan.ovr_rate: + self._ovr_rate[chan.name] = chan.ovr_rate * 100 + else: + self._ovr_rate.pop(chan.name, None) + qual = data.get_quality() + self.input_phase_stddev = qual.input_stddev.imag + self.output_phase_slope = qual.output_slope.imag + self.output_amp_slope = qual.output_slope.real + self.input_amplitude = data.inp.amplitude / 2 ** 15 * self.input_range + self.output_amplitude = data.out.amplitude / 2 ** 15 * self.output_range + self._iq = iq = data.iq * self.output_range / self.input_range 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 + self.amp = np.abs(iq) + self.read_value() + self.set_status(IDLE, '') + elif self._busy: + self._busy = False + if self._action == DONE_GO: + self.set_status(BUSY, 'acquiring') + else: + self.set_status(IDLE, 'acquiring') + return + if self._action == CONTINUE and self.continuous: + self.start_acquisition() + self.set_status(IDLE, 'acquiring') + return + if self._action == GO: + self.start_acquisition() + self._action = DONE_GO + self.set_status(BUSY, 'acquiring') + return + if self._action == DONE_GO: + self._action = WAIT_GO + self._wait_until = time.time() + 2 + self.set_status(IDLE, 'paused') + return + if self._action == WAIT_GO: + if time.time() > self._wait_until: + self._action = CONTINUE + self.start_acquisition() + self.set_status(IDLE, 'acquiring') + def set_status(self, *status): + self._status = status + if self._status != self.status: + self.read_status() + + def read_status(self): + if self._ovr_rate and self._status[0] < WARN: + return WARN, 'overrange on %s' % ' and '.join(self._ovr_rate) + return self._status + + def read_value(self): + if self.imod: + self.imod.value = self._iq.real + if self.qmod: + self.qmod.value = self._iq.imag + return self._iq.real, self._iq.imag + + @Command def go(self): - self.starttime = time.time() - freq = self.freq.value - self._data_args = (RUSdata, freq, self.periods) + """start aquisition""" + if self._busy: + self._action = GO + else: + self._action = DONE_GO + self.start_acquisition() + self._status = BUSY, 'acquiring' + self.read_status() + + def start_acquisition(self): + freq = self.freq.read_value() self.sr = round(self.periods * self.adq.sample_rate / freq) self.adq.init(self.sr, 1) - self.adq.start() - self.read_status() + self.adq.start(RUSdata(self.adq, freq, self.periods)) + self._busy = True + self.setFastPoll(True, 0.001) -class ControlLoop: +class ControlLoop(Module): + roi = Attached(Roi) maxstep = Parameter('max frequency step', FloatRange(unit='Hz'), readonly=False, - default=10000) + 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,