linse/frappy
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

ultrasound: reworked after tests

Browse Source 
- new classes in frappy_psi/ultrasound.py and frappy_psi/adq.mr.py
- add signal plottter
- move clients to bin/ directory

Change-Id: I8db8e5ebc082c346278f09e0e54504e070655f14
This commit is contained in:
zolliker
2025-03-26 15:31:46 +01:00
parent 322cd39e0a
commit 6f547f0781
7 changed files with 553 additions and 167 deletions
Download Patch File Download Diff File Expand all files Collapse all files

110
frappy_psi/adq_mr.py
View File

@ -50,6 +50,27 @@ GHz = 1e9
RMS_TO_VPP = 2 * np.sqrt(2)
class Timer:
def __init__(self):
self.data = [(time.time(), 'start')]
def __call__(self, text=''):
now = time.time()
prev = self.data[-1][0]
self.data.append((now, text))
return now - prev
def summary(self):
return ' '.join(f'{txt} {tim:.3f}' for tim, txt in self.data[1:])
def show(self):
first = prev = self.data[0][0]
print('---', first)
for tim, txt in self.data[1:]:
print(f'{(tim - first) * 1000:9.3f} {(tim - prev) * 1000:9.3f} ms {txt}')
prev = tim
class Adq:
sample_rate = 2 * GHz
max_number_of_channels = 2
@ -161,7 +182,6 @@ class Adq:
ADQAPI.ADQ_DisarmTrigger(self.adq_cu, self.adq_num)
ADQAPI.ADQ_ArmTrigger(self.adq_cu, self.adq_num)
self.data = data
self.starttime = time.time()
def get_data(self):
"""get new data if available"""
@ -173,12 +193,14 @@ class Adq:
if ADQAPI.ADQ_GetAcquiredAll(self.adq_cu, self.adq_num):
ready = True
data.timer('ready')
else:
if self.trigger == SW_TRIG:
ADQAPI.ADQ_SWTrig(self.adq_cu, self.adq_num)
if not ready:
self.busy = True
return None
self.data = None
t = time.time()
# Get data from ADQ
if not ADQAPI.ADQ_GetData(
@ -187,7 +209,6 @@ class Adq:
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
@ -197,12 +218,20 @@ class PEdata:
self.sample_rate = adq.sample_rate
self.samp_freq = self.sample_rate / GHz
self.number_of_records = adq.number_of_records
self.timer = Timer()
def retrieve(self, adq):
data = []
rawsignal = []
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
rawsignal.append(onedim[:adq.samples_per_record])
# convert 16 bit int to a value in the range -1 .. 1
data.append(onedim.reshape(adq.number_of_records, adq.samples_per_record) / float(2 ** 15))
# Now this is an array with all records, but the time is artificial
self.data = data
self.rawsignal = rawsignal
self.timer('retrieved')
def sinW(self, sig, freq, ti, tf):
# sig: signal array
@ -232,6 +261,7 @@ 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)])
@ -254,24 +284,32 @@ class PEdata:
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()))
self.timer('gates')
try:
self.ndecimate = int(round(self.sample_rate / freq))
except TypeError as e:
raise TypeError(f'{self.sample_rate}/{freq} {e}')
iq = self.averageiq(self.data, freq / GHz, *pulse)
# times.append(('filtro', time.time()))
self.timer('aviq')
iqf = self.filtro(iq, bw_cutoff)
m = len(iqf[0]) // self.ndecimate
self.timer('filtro')
m = max(1, len(iqf[0]) // self.ndecimate)
ll = m * self.ndecimate
iqf = [iqfx[0:ll] for iqfx in iqf]
# times.append(('iqdec', time.time()))
self.timer('iqf')
iqd = np.average(np.resize(iqf, (2, m, self.ndecimate)), axis=2)
self.timer('avg')
t_axis = np.arange(m) * self.ndecimate / self.samp_freq
pulsig = np.abs(self.data[0][0])
# times.append(('pulsig', time.time()))
self.timer('pulsig')
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]
result = ([self.box(iqf, *r) for r in roi], # gates
(t_axis, iqd[0], iqd[1], pulsig)) # curves
self.timer('result')
# self.timer.show()
# ns = len(self.rawsignal[0]) * self.number_of_records
# print(f'{ns} {ns / 2e6} ms')
return result
class Namespace:
@ -290,24 +328,40 @@ class RUSdata:
self.inp = Namespace(idx=0, name='input')
self.out = Namespace(idx=1, name='output')
self.channels = (self.inp, self.out)
self.timer = Timer()
def retrieve(self, adq):
self.timer('start retrieve')
npts = self.samples_per_record - self.delay_samples
complex_sinusoid = np.exp(1j * 2 * np.pi * self.freq / self.sample_rate * np.arange(npts))
nbin = max(1, npts // (self.periods * 60)) # for performance reasons, do the binning first
nreduced = npts // nbin
ft = 2 * np.pi * self.freq * nbin / self.sample_rate * np.arange(nreduced)
self.timer('create time axis')
# complex_sinusoid = np.exp(1j * ft) # do not use this, below is 33 % faster
complex_sinusoid = 1j * np.sin(ft) + np.cos(ft)
self.timer('sinusoid')
rawsignal = [] # for raw plot
for chan in self.channels: # looping over input and output
# although the ADC is only 14 bit it is representet as unsigend 16 bit numbers,
# although the ADC is only 14 bit it is represented as unsigend 16 bit numbers,
# and due to some calculations (calibration) the last 2 bits are not zero
isignal = np.frombuffer(adq.target_buffers[chan.idx].contents, dtype=np.int16)[self.delay_samples:]
# importatn convert to float first, before doing any calculations.
# calculations with int16 may silently overflow
chan.signal = isignal / float(2 ** 16) # in V -> peak to peak 1 V ~ +- 0.5 V
chan.ovr_rate = (np.abs(isignal) > 32000).sum() / npts # fraction of points near end of range
beg = self.delay_samples
isignal = np.frombuffer(adq.target_buffers[chan.idx].contents, dtype=np.int16)[beg:beg+nreduced * nbin]
self.timer('isignal')
reduced = isignal.reshape((-1, nbin)).mean(axis=1) # this converts also int16 to float
self.timer('reduce')
rawsignal.append(reduced)
chan.signal = signal = reduced * 2 ** -16 # in V -> peak to peak 1 V ~ +- 0.5 V
self.timer('divide')
# calculate RMS * sqrt(2) -> peak sinus amplitude.
# may be higher than the input range by a factor 1.4 when heavily clipped
signal = chan.signal - np.median(chan.signal)
chan.amplitude = np.sqrt((signal ** 2).mean()) * RMS_TO_VPP
self.timer('amp')
chan.mixed = signal * complex_sinusoid
self.timer('mix')
chan.mean = chan.mixed.mean()
self.timer('mean')
self.rawsignal = rawsignal
if self.inp.mean:
self.iq = self.out.mean / self.inp.mean
else:
@ -324,15 +378,19 @@ class RUSdata:
the imaginary part indicates a turning phase (rad/sec)
the real part indicates changes in amplitude (0.01 ~= 1%/sec)
"""
self.timer('get_quality')
npts = len(self.channels[0].signal)
nper = npts // self.periods
nbin = 50 # 50 samples, 25 ns
for chan in self.channels:
mean = chan.mixed.mean()
chan.reduced = chan.mixed[:self.periods * nper].reshape((-1, nper)).mean(axis=1) / mean
chan.binned = chan.signal[:npts // nbin * nbin].reshape((-1, nbin)).mean(axis=1)
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])
result = Namespace(
input_stddev=self.inp.reduced.std(),
output_slope=np.polyfit(timeaxis, self.out.reduced, 1)[0])
self.timer('got_quality')
self.timer.show()
ns = len(self.rawsignal[0])
print(f'{ns} {ns / 2e6} ms')
return result