From 2e99e45aeaef6cacbb35a7c97f5e73e7eb7b09a5 Mon Sep 17 00:00:00 2001 From: Markus Zolliker Date: Tue, 28 Jan 2025 09:40:36 +0100 Subject: [PATCH] WIP new version of ultrasound Change-Id: Iadb83396a64e277f6f0a37f7a96d92105648c4fe --- frappy_psi/adq_mr.py | 569 +++++++++++++++++++++------------------ frappy_psi/ultrasound.py | 374 +++++++++++++------------ 2 files changed, 517 insertions(+), 426 deletions(-) diff --git a/frappy_psi/adq_mr.py b/frappy_psi/adq_mr.py index 1f85b11..787ba9c 100644 --- a/frappy_psi/adq_mr.py +++ b/frappy_psi/adq_mr.py @@ -1,256 +1,313 @@ -""" -Created on Tue Nov 26 15:42:43 2019 - -@author: tartarotti_d-adm -""" - - -import sys -import atexit -import signal -import time -import numpy as np -import ctypes as ct -from numpy import sqrt, arctan2, sin, cos -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 -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 = scipy.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() - 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 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""" - #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 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(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 - - 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""" - 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 - 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] - +# ***************************************************************************** +# 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 +# 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 +# ***************************************************************************** +"""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 scipy.signal import butter, filtfilt + + +# For different trigger modes +SW_TRIG = 1 +# 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 +INT_TRIG = 4 +LVL_FALLING = 0 +LVL_RISING = 1 + +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 + +GHz = 1e9 + + +class Adq: + sample_rate = 2 * GHz + max_number_of_channels = 2 + ndecimate = 50 # decimation ratio (2GHz / 40 MHz) + number_of_records = 1 + 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") + + 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, '.') + + # Find ADQ devices + 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) + + 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('') + + 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 trigger mode + if not ADQAPI.ADQ_SetTriggerMode(self.adq_cu, self.adq_num, self.trigger): + raise RuntimeError('ADQ_SetTriggerMode failed.') + if self.trigger == LVL_TRIG: + if not ADQAPI.ADQ_SetLvlTrigLevel(self.adq_cu, self.adq_num, -100): + raise RuntimeError('ADQ_SetLvlTrigLevel failed.') + if not ADQAPI.ADQ_SetTrigLevelResetValue(self.adq_cu, self.adq_num, 1000): + raise RuntimeError('ADQ_SetTrigLevelResetValue failed.') + if not ADQAPI.ADQ_SetLvlTrigChannel(self.adq_cu, self.adq_num, 1): + raise RuntimeError('ADQ_SetLvlTrigChannel failed.') + if not ADQAPI.ADQ_SetLvlTrigEdge(self.adq_cu, self.adq_num, LVL_RISING): + raise RuntimeError('ADQ_SetLvlTrigEdge failed.') + elif self.trigger == EXT_TRIG_1: + if not ADQAPI.ADQ_SetExternTrigEdge(self.adq_cu, self.adq_num, 2): + raise RuntimeError('ADQ_SetLvlTrigEdge failed.') + # 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""" + 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)() + 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 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) + self.status = self.BUSY + + 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, 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 + + +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(adq.target_buffers[ch].contents, dtype=np.int16) + data.append(onedim.reshape(adq.number_of_records, adq.samples_per_record) / float(2**14)) # 14 bits ADC + # Now this is an array with all records, but the time is artificial + self.data = data + + 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 of sigin + a, b = self.sinW(sigin, freq, ti, tf) + amp = np.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 * 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)]) + return iorq.sum(axis=0) / self.number_of_records + + def filtro(self, iorq, cutoff): + # butter lowpass + nyq = 0.5 * self.sample_rate + normal_cutoff = cutoff / nyq + order = 5 + 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): + 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, freq, pulse, roi, bw_cutoff): + """return iq values of rois and prepare plottable curves for iq""" + self.ndecimate = int(round(self.sample_rate / freq)) + # times = [] + # times.append(('aviq', time.time())) + iq = self.averageiq(self.data, freq / GHz, *pulse) + # times.append(('filtro', time.time())) + 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())) + 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(self.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] + + +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] diff --git a/frappy_psi/ultrasound.py b/frappy_psi/ultrasound.py index a35559e..ec46528 100644 --- a/frappy_psi/ultrasound.py +++ b/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 +from frappy_psi.adq_mr import Adq, PEdata, RUSdata from frappy.core import Attached, BoolType, Done, FloatRange, HasIO, \ - IntRange, Module, Parameter, Readable, StringIO, StringType, \ - IDLE, DISABLED, TupleOf, ArrayOf + IntRange, Module, Parameter, Readable, Writable, Drivable, StringIO, StringType, \ + 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): @@ -55,7 +56,7 @@ class Roi(Readable): enable = Parameter('calculate this roi', BoolType(), readonly=False, default=True) pollinterval = Parameter(export=False) - interval = (0,0) + interval = (0, 0) def initModule(self): super().initModule() @@ -92,80 +93,22 @@ class FreqStringIO(StringIO): end_of_line = '\r' -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) +class Frequency(HasIO, Writable): + value = Parameter('frequency', unit='Hz') 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 - old = None - starttime = None - interval = (0,0) + def register_dif(self, dif): + self.dif = dif - def earlyInit(self): - super().earlyInit() - 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() - - 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_target(self, value): + self.communicate('FREQ %.15g;FREQ?' % value) + self.last_change = time.time() + if self.dif: + self.dif.read_value() 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): - """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() +class FrequencyDif(Readable): + freq = Attached(Frequency) + base = Parameter('base frequency', FloatRange(unit='Hz'), default=0) + value = Parameter('difference to base frequency', FloatRange(unit='Hz'), default=0) - if self.starttime is None: - self.starttime = time.time() - 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] + def initModule(self): + super().initModule() + self.freq.register_dif(self) - 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 - 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 - if self.control: - newfreq = freq + inphase * self.slope - self.basefreq - # step = sorted((-self.maxstep, inphase * self.slope, self.maxstep))[1] - if self.old: - fdif = freq - self.old[0] - idif = inphase - self.old[1] - if abs(fdif) >= self.minstep: - self.slope = - fdif / idif - else: - fdif = 0 - idif = 0 - newfreq = freq + self.minstep - self.old = (freq, inphase) - if self.skipctrl > 0: # do no control for some time after changing frequency - self.skipctrl -= 1 - elif self.control: - self.freq = sorted((self.freq - self.maxstep, newfreq, self.freq + self.maxstep))[1] - #print(times) - return Done + def read_value(self): + return self.freq - self.base -class Curves(Readable): +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() + + def read_value(self): + if self.get_rawdata(): # new data available + roilist = [r for r in self.roilist if r.enable] + freq = self.freq.value + gates = self.adq.gates_and_curves(self._data, freq, + (self.time, self.time + self.size), + [r.interval for r in roilist]) + 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 + return self.adq.curves + + # TODO: CONTROL + # inphase = self.roilist[0].i + # if self.control: + # newfreq = freq + inphase * self.slope - self.basefreq + # # step = sorted((-self.maxstep, inphase * self.slope, self.maxstep))[1] + # if self.old: + # fdif = freq - self.old[0] + # idif = inphase - self.old[1] + # if abs(fdif) >= self.minstep: + # self.slope = - fdif / idif + # else: + # fdif = 0 + # idif = 0 + # newfreq = freq + self.minstep + # self.old = (freq, inphase) + # if self.skipctrl > 0: # do no control for some time after changing frequency + # self.skipctrl -= 1 + # elif self.control: + # 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()) + + 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))