Source code for damaris.data.ADC_Result

# -*- coding: iso-8859-1 -*-

from __future__ import annotations

from .Resultable import Resultable
from .Drawable import Drawable
from .Signalpath import Signalpath
from .DamarisFFT import DamarisFFT
import threading
import numpy
import sys
import datetime
import tables
from typing import IO
from scipy.signal import filtfilt, remez

_XArray = numpy.ndarray | list
_YChannels = list[numpy.ndarray]
_IndexList = list[tuple[int, int]]
_Description = dict[str, str]
#############################################################################
#                                                                           #
# Name: Class ADC_Result                                                    #
#                                                                           #
# Purpose: Specialised class of Resultable and Drawable                     #
#          Contains recorded ADC Data                                       #
#                                                                           #
#############################################################################

[docs] class ADC_Result(Resultable, Drawable, DamarisFFT, Signalpath): default_bit_depth = 14 """ Represents the result of an ADC, encapsulating data and metadata for processing, visualization, and export. This class combines data storage and manipulation functionality with interfaces for drawing and result processing. It manages time-series data across multiple channels, supports dynamic resizing of datasets, and provides mechanisms to export data in various formats. Its key roles include ADC result storage, metadata management, and integration with external systems. Attributes: xlabel: Label for the x-axis, used in visualization. ylabel: Label for the y-axis, used in visualization. lock: A threading lock for synchronizing access to the data. nChannels: Number of data channels in the result set. sampling_rate: Sampling frequency of the ADC data. job_id: identifier for the job that generated this result. """
[docs] def __init__(self, x: _XArray | None = None, y: _YChannels | None = None, index: _IndexList | None = None, sampl_freq: float | None = None, desc: _Description | None = None, job_id: int | None = None, job_date: datetime.datetime | None = None): Resultable.__init__(self) Drawable.__init__(self) # Title of this accumulation: set Values: Job-ID and Description (plotted in GUI -> look Drawable) # Is set in ResultReader.py (or in copy-construktor) self.__title_pattern = "ADC-Result: job_id = %s, desc = %s" # Axis-Labels (inherited from Drawable) self.xlabel = "Time (s)" self.ylabel = "Samples [Digits]" self.lock=threading.RLock() self.nChannels = 0 self.is_clipped = False # using no argument for initialization if (x is None) and (y is None) and (index is None) and (sampl_freq is None) and (desc is None) and (job_id is None) and (job_date is None): self.cont_data = False self.sampling_rate = 0 self.index = [] self.x = numpy.array([], dtype="float32") self.y = [] # using all arguments for initialization elif (x is not None) and (y is not None) and (index is not None) and (sampl_freq is not None) and (desc is not None) and (job_id is not None) and (job_date is not None): # TODO: insure integer calculations for ADC_Result operations. #for ch in y: # if not numpy.issubdtype(ch.dtype, numpy.integer): # raise TypeError("TypeError: ADC_Result y data must be a list with integer type channels") self.x = numpy.asarray(x) self.y = [numpy.asarray(ch) for ch in y] self.index = index self.sampling_rate = sampl_freq self.cont_data = True self.description = desc self.job_id = job_id self.job_date = job_date title="ADC-Result: job-id=%d"%int(self.job_id) if len(self.description)>0: for k,v in self.description.items(): # string keys can be made invisible by adding two underscores in front of them if not (type(k) in (str,) and k[0] == '_' and k[1] == '_'): title+=", %s=%s"%(k,v) self.set_title(title) else: raise ValueError("Wrong usage of __init__!")
def create_data_space(self, channels: int, samples: int) -> None: "Initialises the internal data-structures" if self.contains_data(): print("Warning ADC-Result: Tried to run \"create_data_space()\" more than once.") return for i in range(channels): self.y.append(numpy.zeros((samples,), dtype="int16")) self.x = numpy.zeros((samples,), dtype="float32") self.nChannels = channels self.index.append((0, samples-1)) self.cont_data = True def contains_data(self) -> bool: """Returns true if ADC_Result contains data. (-> create_data_space() was called)""" return self.cont_data def check_clipping(self, bit_depth: int | None = None) -> bool: """ Check if any data in y is clipped. For a 14-bit ADC, clipping occurs at 2**13 (8192) or 2**13-1 (8191). In general, for a bit_depth-bit signed ADC, limits are 2**(bit_depth-1) and -(2**(bit_depth-1)). """ if not self.contains_data(): return False if bit_depth is None: # try to get bit_depth from description if hasattr(self, "description") and self.description is not None: if "adc_bit_depth" in self.description: try: bit_depth = int(self.description["adc_bit_depth"]) except (ValueError, TypeError): bit_depth = ADC_Result.default_bit_depth else: bit_depth = ADC_Result.default_bit_depth else: bit_depth = ADC_Result.default_bit_depth limit_pos = 2**(bit_depth-1) - 1 limit_neg = -2**(bit_depth-1) self.lock.acquire() try: self.is_clipped = False for ch_data in self.y: if numpy.any(ch_data >= limit_pos) or numpy.any(ch_data <= limit_neg): self.is_clipped = True break finally: self.lock.release() return self.is_clipped def lowpass(self, cutoff: float, numtaps: int | None = None) -> ADC_Result: """ Apply a zero-phase lowpass FIR filter to all channels. Uses scipy.signal.filtfilt (forward-backward filtering) with a Parks-McClellan FIR filter (remez). The forward-backward pass guarantees zero phase distortion and doubles the effective filter order. The stopband is automatically set to 0.1 * Nyquist beyond the cutoff. If numtaps is not provided, it is estimated from the transition width to yield approximately 40-50 dB stopband attenuation (typically ~40 taps). Parameters ---------- cutoff : float Passband edge frequency in Hz. Frequencies above this value are attenuated. Must be a positive number below the Nyquist frequency. numtaps : int, optional Number of FIR taps. If omitted, computed automatically from the transition width. Returns ------- ADC_Result A new ADC_Result with the filtered data (float64 per channel). The original object is not modified. Raises ------ RuntimeError If the result contains no data. ValueError If cutoff is invalid for the current sampling rate. Examples -------- >>> # Lowpass at 1 MHz (taps auto-computed, ~40) ... filtered = adc.lowpass(cutoff=1e6) >>> # Lowpass at 500 kHz with fixed 63 taps ... filtered = adc.lowpass(cutoff=500e3, numtaps=63) """ if not self.contains_data(): raise RuntimeError("lowpass: no data present") nyquist = self.sampling_rate / 2.0 cutoff = float(cutoff) if cutoff <= 0 or cutoff >= nyquist: raise ValueError( f"cutoff {cutoff} Hz is outside valid range (0, {nyquist})" ) # Transition width: 0.1 * Nyquist trans_width = 0.1 * nyquist stopband = min(cutoff + trans_width, nyquist * 0.99) # Auto-estimate taps if not given if numtaps is None: numtaps = max(21, int(round(8.0 / (stopband - cutoff) * nyquist))) if numtaps % 2 == 0: numtaps += 1 coeffs = remez(numtaps, [0, cutoff, stopband, nyquist], [1, 0], fs=self.sampling_rate) self.lock.acquire() try: tmp_y = [] for i in range(self.get_number_of_channels()): data = numpy.asarray(self.y[i], dtype='float64') filtered = filtfilt( coeffs, [1.0], data, padtype='odd', padlen=(numtaps - 1) * 3 ) tmp_y.append(filtered) r = ADC_Result( x=self.x[:], y=tmp_y, index=self.index[:], sampl_freq=self.sampling_rate, desc=self.description, job_id=self.job_id, job_date=self.job_date, ) self._copy_metadata(r) finally: self.lock.release() return r def add_sample_space(self, samples: int) -> None: "Adds space for n samples, where n can also be negative (deletes space). New space is filled up with \"0\"" self.lock.acquire() if not self.cont_data: print("Warning ADC-Result: Tried to resize empty array!") return length = len(self.y[0]) self.x = numpy.resize(self.x, (length+samples)) for i in range(self.get_number_of_channels()): self.y[i] = numpy.resize(self.y[i], (length+samples)) self.index.append((length, len(self.y[0])-1)) self.lock.release() def get_result_by_index(self, index: int) -> ADC_Result: self.lock.acquire() try: start = self.index[index][0] end = self.index[index][1] except: self.lock.release() raise tmp_x = self.x[start:end+1].copy() tmp_y = [] for i in range(self.get_number_of_channels()): tmp_y.append(self.y[i][start:end+1].copy()) r = ADC_Result(x = tmp_x, y = tmp_y, index = [(0,len(tmp_y[0])-1)], sampl_freq = self.sampling_rate, desc = self.description.copy() if self.description is not None else None, job_id = self.job_id, job_date = self.job_date) self._copy_metadata(r) self.lock.release() return r def get_sampling_rate(self) -> float: """Returns the samplingfrequency""" return self.sampling_rate + 0 def set_sampling_rate(self, hz: float) -> None: """Sets the samplingfrequency in hz""" self.sampling_rate = float(hz) def get_nChannels(self) -> int: """Gets the number of channels""" return self.nChannels + 0 def set_nChannels(self, channels: int) -> None: """Sets the number of channels""" self.nChannels = int(channels) def get_index_bounds(self, index: int) -> tuple[int, int]: "Returns a tuple with (start, end) of the wanted result" return self.index[index] def uses_statistics(self) -> bool: return False def write_to_csv(self, destination: str | IO[str] = sys.stdout, delimiter: str = " ") -> None: """ writes the data to a file or to sys.stdout destination can be a file or a filename suitable for further processing """ if isinstance(destination, str): with open(destination, "w") as f: self._write_to_csv_dest(f, delimiter) else: self._write_to_csv_dest(destination, delimiter) def _write_to_csv_dest(self, dest: IO[str], delimiter: str) -> None: dest.write("# adc_result\n") dest.write("# t y0 y1 ...\n") self.lock.acquire() try: xdata=self.get_xdata() ch_no=self.get_number_of_channels() ydata=list(map(self.get_ydata, range(ch_no))) for i in range(len(xdata)): dest.write("%e"%xdata[i]) for j in range(ch_no): dest.write("%s%e"%(delimiter, ydata[j][i])) dest.write("\n") finally: self.lock.release() def write_to_simpson(self, destination: str | IO[str] = sys.stdout, delimiter: str = " ") -> None: """ writes the data to a text file or sys.stdout in Simpson format, for further processing with the NMRnotebook software; destination can be a file or a filename """ if isinstance(destination, str): with open(destination, "w") as f: self._write_to_simpson_dest(f, delimiter) else: self._write_to_simpson_dest(destination, delimiter) def _write_to_simpson_dest(self, dest: IO[str], delimiter: str) -> None: self.lock.acquire() try: xdata=self.get_xdata() dest.write("SIMP\n") dest.write("%s%i%s"%("NP=", len(xdata), "\n")) dest.write("%s%i%s"%("SW=", self.get_sampling_rate(), "\n")) dest.write("TYPE=FID\n") dest.write("DATA\n") ch_no=self.get_number_of_channels() ydata=list(map(self.get_ydata, range(ch_no))) for i in range(len(xdata)): for j in range(ch_no): dest.write("%g%s"%(ydata[j][i], delimiter)) dest.write("\n") dest.write("END\n") finally: self.lock.release() def write_to_hdf(self, hdffile: tables.File, where: str, name: str, title: str, complib: str | None = None, complevel: int | None = None) -> None: accu_group=hdffile.create_group(where=where,name=name,title=title) accu_group._v_attrs.damaris_type="ADC_Result" if self.contains_data(): self.lock.acquire() try: # save time stamps if "job_date" in dir(self) and self.job_date is not None: accu_group._v_attrs.time="%04d%02d%02d %02d:%02d:%02d.%03d"%(self.job_date.year, self.job_date.month, self.job_date.day, self.job_date.hour, self.job_date.minute, self.job_date.second, self.job_date.microsecond/1000) if self.description is not None: for (key,value) in self.description.items(): accu_group._v_attrs.__setattr__("description_"+key,str(value)) accu_group._v_attrs.__setattr__("sampling_rate",self.sampling_rate) # save interval information filter=None if complib is not None: if complevel is None: complevel=9 filter=tables.Filters(complevel=complevel,complib=complib,shuffle=1) index_table=hdffile.create_table(where=accu_group, name="indices", description={"start": tables.UInt64Col(), "length": tables.UInt64Col(), "start_time": tables.Float32Col(), "dwelltime": tables.Float32Col()}, title="indices of adc data intervals", filters=filter, expectedrows=len(self.index)) index_table.flavor="numpy" # save channel data new_row=index_table.row for i in range(len(self.index)): new_row["start"]=self.index[i][0] new_row["dwelltime"]=1.0/self.sampling_rate new_row["start_time"]=1.0/self.sampling_rate*self.index[i][0] new_row["length"]=self.index[i][1]-self.index[i][0]+1 new_row.append() index_table.flush() new_row=None index_table=None # prepare saving data channel_no=len(self.y) timedata=numpy.empty((len(self.y[0]),channel_no), dtype = "int16") for ch in range(channel_no): timedata[:,ch]=self.get_ydata(ch) # save data time_slice_data=None if filter is not None: chunkshape = numpy.shape(timedata) if len(chunkshape) <= 1: chunkshape = (min(chunkshape[0],1024*8),) else: chunkshape = (min(chunkshape[0],1024*8), chunkshape[1]) if tables.__version__[0]=="1": time_slice_data=hdffile.create_carray(accu_group, name="adc_data", shape=timedata.shape, atom=tables.Int16Atom(shape=chunkshape, flavor="numpy"), filters=filter, title="adc data") else: time_slice_data=hdffile.create_carray(accu_group, name="adc_data", shape=timedata.shape, chunkshape=chunkshape, atom=tables.Int16Atom(), filters=filter, title="adc data") time_slice_data[:]=timedata else: time_slice_data=hdffile.create_array(accu_group, name="adc_data", obj=timedata, title="adc data") finally: timedata=None time_slice_data=None accu_group=None self.lock.release() # Ueberladen von Operatoren und Built-Ins ------------------------------------------------------- def __len__(self) -> int: "Redefining len(ADC_Result obj), returns the number of samples in one channel and 0 without data" if len(self.y)>0: return len(self.y[0]) return 0 def __repr__(self) -> str: """ writes job meta data and data to string returned """ tmp_string = "Job ID: " + str(self.job_id) + "\n" tmp_string += "Job Date: " + str(self.job_date) + "\n" tmp_string += "Description: " + str(self.description) + "\n" if len(self.y)>0: tmp_string += "Indexes: " + str(self.index) + "\n" tmp_string += "Samples per Channel: " + str(len(self.y[0])) + "\n" tmp_string += "Samplingfrequency: " + str(self.sampling_rate) + "\n" tmp_string += "X: " + repr(self.x) + "\n" for i in range(self.get_number_of_channels()): tmp_string += ("Y(%d): " % i) + repr(self.y[i]) + "\n" return tmp_string def _copy_metadata(self, r: ADC_Result) -> None: "Copy channel metadata (nChannels, is_clipped) from self to the new result r" r.nChannels = self.nChannels r.is_clipped = self.is_clipped def __add__(self, other: int | float) -> ADC_Result: "Redefining self + other (scalar)" if isinstance(other, int) or isinstance(other, float): self.lock.acquire() tmp_y = [] for i in range(self.get_number_of_channels()): tmp_y.append(numpy.array(self.y[i], dtype="float32") + other) r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date) self._copy_metadata(r) self.lock.release() return r else: raise ValueError(f"ValueError: Cannot add \"{other.__class__}\" to ADC-Result!") def __radd__(self, other: int | float) -> ADC_Result: "Redefining other (scalar) + self" return self.__add__(other) def __sub__(self, other: int | float) -> ADC_Result: "Redefining self - other (scalar)" if isinstance(other, int) or isinstance(other, float): self.lock.acquire() tmp_y = [] for i in range(self.get_number_of_channels()): tmp_y.append(numpy.array(self.y[i], dtype="float32") - other) r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date) self._copy_metadata(r) self.lock.release() return r else: raise ValueError(f"ValueError: Cannot subtract \"{other.__class__}\" to ADC-Result!") def __rsub__(self, other: int | float) -> ADC_Result: "Redefining other (scalar) - self" if isinstance(other, int) or isinstance(other, float): self.lock.acquire() tmp_y = [] for i in range(self.get_number_of_channels()): tmp_y.append(other - numpy.array(self.y[i], dtype="float32")) r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date) self._copy_metadata(r) self.lock.release() return r else: raise ValueError(f"ValueError: Cannot subtract \"{other.__class__}\" to ADC-Result!") def __mul__(self, other: int | float) -> ADC_Result: "Redefining self * other (scalar)" if isinstance(other, int) or isinstance(other, float): self.lock.acquire() tmp_y = [] for i in range(self.get_number_of_channels()): tmp_y.append(numpy.array(self.y[i], dtype="float32") * other) r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date) self._copy_metadata(r) self.lock.release() return r else: raise ValueError(f"ValueError: Cannot multiply \"{other.__class__}\" to ADC-Result!") def __rmul__(self, other: int | float) -> ADC_Result: "Redefining other (scalar) * self" return self.__mul__(other) def __pow__(self, other: int | float) -> ADC_Result: "Redefining self ** other (scalar)" if isinstance(other, int) or isinstance(other, float): self.lock.acquire() tmp_y = [] for i in range(self.get_number_of_channels()): tmp_y.append(numpy.array(self.y[i], dtype="float32") ** other) r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date) self._copy_metadata(r) self.lock.release() return r else: raise ValueError(f"ValueError: Cannot power raise \"{other.__class__}\" to ADC-Result!") def __truediv__(self, other: int | float) -> ADC_Result: "Redefining self / other (scalar)" if isinstance(other, int) or isinstance(other, float): self.lock.acquire() tmp_y = [] for i in range(self.get_number_of_channels()): tmp_y.append(numpy.array(self.y[i], dtype="float32") / other) r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date) self._copy_metadata(r) self.lock.release() return r else: raise ValueError(f"ValueError: Cannot divide \"{other.__class__}\" to ADC-Result!") def __rtruediv__(self, other: int | float) -> ADC_Result: "Redefining other (scalar) / self" if isinstance(other, int) or isinstance(other, float): self.lock.acquire() tmp_y = [] for i in range(self.get_number_of_channels()): tmp_y.append(other / numpy.array(self.y[i], dtype="float32")) r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date) self._copy_metadata(r) self.lock.release() return r else: raise ValueError(f"ValueError: Cannot divide \"{other.__class__}\" to ADC-Result!") def __floordiv__(self, other: int) -> ADC_Result: "Redefining self // other (scalar)" if isinstance(other, float): raise ValueError("ValueError: Cannot use floor division (//) on floats! Use \"/\" instead of \"//\"! ") if isinstance(other, int): self.lock.acquire() tmp_y: list[numpy.ndarray] = [] for i in range(self.get_number_of_channels()): tmp_y.append(numpy.array(self.y[i], dtype="float32") // other) r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date) self._copy_metadata(r) self.lock.release() return r else: raise ValueError(f"ValueError: Cannot divide \"{other.__class__}\" to ADC-Result!") def __rfloordiv__(self, other: int) -> ADC_Result: "Redefining other (scalar) // self" if isinstance(other, float): raise ValueError("ValueError: Cannot use floor division (//) on floats! Use \"/\" instead of \"//\"! ") if isinstance(other, int): self.lock.acquire() tmp_y: list[numpy.ndarray] = [] for i in range(self.get_number_of_channels()): tmp_y.append(other // numpy.array(self.y[i], dtype="float32")) r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date) self._copy_metadata(r) self.lock.release() return r else: raise ValueError(f"ValueError: Cannot divide \"{other.__class__}\" to ADC-Result!") def __neg__(self) -> ADC_Result: "Redefining -self" self.lock.acquire() tmp_y = [] for i in range(self.get_number_of_channels()): tmp_y.append(numpy.array(-self.y[i])) r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date) self._copy_metadata(r) self.lock.release() return r
def read_from_hdf(hdf_node: object) -> ADC_Result | None: """ read accumulation data from HDF node and return it. """ # formal checks first if not isinstance(hdf_node, tables.Group): return None if hdf_node._v_attrs.damaris_type!="ADC_Result": return None if not (hdf_node.__contains__("indices") and hdf_node.__contains__("adc_data")): return None # job id and x,y titles are missing adc=ADC_Result() if hdf_node._v_title: adc.set_title(hdf_node._v_title) # populate description dictionary adc.description={} for attrname in hdf_node._v_attrs._v_attrnamesuser: if attrname.startswith("description_"): adc.description[attrname[12:]]=hdf_node._v_attrs.__getattr__(attrname) if "time" in dir(hdf_node._v_attrs): timestring=hdf_node._v_attrs.__getattr__("time") adc.job_date=datetime.datetime(int(timestring[:4]), # year int(timestring[4:6]), # month int(timestring[6:8]), # day int(timestring[9:11]), # hour int(timestring[12:14]), # minute int(timestring[15:17]), # second int(timestring[18:21])*1000 # microsecond ) # start with indices for r in hdf_node.indices.iterrows(): adc.index.append((r["start"],r["start"]+r["length"]-1)) adc.sampling_rate=1.0/r["dwelltime"] # now really belief there are no data if len(adc.index)==0: adc.cont_data=False return adc adc.cont_data=True # now do the real data adc_data=hdf_node.adc_data.read() adc.set_nChannels(adc_data.shape[1]) adc.x=numpy.arange(adc_data.shape[0], dtype="float32")/adc.sampling_rate for ch in range(adc_data.shape[1]): adc.y.append(adc_data[:,ch]) adc.check_clipping() return adc