Patch13 dev

docsrc/API_docs.rst

@@ -36,6 +36,7 @@ Sequences
    pyepr.sequences.CarrPurcellSequence
    pyepr.sequences.ResonatorProfileSequence
    pyepr.sequences.TWTProfileSequence
+   pyepr.sequences.T1InversionRecoverySequence
 
 Pulses
 ~~~~~~
@@ -90,6 +91,7 @@ I/O
     pyepr.dataset.create_dataset_from_sequence
     pyepr.dataset.create_dataset_from_axes
     pyepr.dataset.create_dataset_from_bruker
+    pyepr.dataset.downconvert_dataset
 
 Utilities
 ~~~~~~~~~

docsrc/install.rst

@@ -38,5 +38,4 @@ PyEPR requires:
     - h5netcdf
     - toml
     - deerlab (https://github.com/JeschkeLab/DeerLab)
-    - numba
     - psutil

docsrc/releasenotes.rst

@@ -1,6 +1,19 @@
 Release Notes
 =============
 
+Version 1.1 (2026-05-27):
+++++++++++++++++++++++++++++
+- Added `AmplifierLinearityAnalysis` class for characterizing amplifier non-linearity.
+- Added `T1InversionRecovery` sequence.
+- Added rise time to linear chirp pulses.
+- Added right based arithmatic.
+- Fixed Version Detection Bug
+- Fixed Dependency issues
+- Improved Documentation
+- Removed Numba dependency
+
+
+
 Version 1.0.0 (2025-09-12):
 ++++++++++++++++++++++++++++
 - All references to `LO` have been changed to `freq` in the frequency object and related.

docsrc/tutorial_sequencer.md

@@ -1,6 +1,6 @@
 # Pulse Sequencer
 
-PyEPR provides an intuitive object-oriented pulse programmer allowing the user to design pulsesequences in a hardware-agnostic manner. Additionally, several common EPR experiments are pre-defined and can be easily instantiated and modified.
+PyEPR provides an intuitive object-oriented pulse programmer allowing the user to design pulse sequences in a hardware-agnostic manner. Additionally, several common EPR experiments are pre-defined and can be easily instantiated and modified.
 
 PyEPR uses ns, GHz and G as the default time, frequency and field units. Very occasionally, other units such as µs or MHz are used, in which case it will be explicitly mentioned.
 
@@ -31,33 +31,33 @@ These pulses will eventually need a scale (amplitude), before the sequence can b
 A Detection window is also created
 ```python
 p90 = epr.RectPulse(tp=16,
-                    freq=0, # Frequency offset in MHz, w.r.t the sequence frequency,
+                    freq=0, # Frequency offset in GHz, w.r.t the sequence frequency,
                     flipangle=np.pi/2, # Flip angle in degrees
                     pcyc = {"phases":[0, np.pi], "dets":[1,-1]}
                     )
 p180 = epr.RectPulse(tp=32,
-                    freq=0, # Frequency offset in MHz, w.r.t the sequence frequency,
+                    freq=0, # Frequency offset in GHz, w.r.t the sequence frequency,
                     flipangle=np.pi, # Flip angle in degrees
                     )
-det = epr.Detetction(tp=32,
-                    freq=0, # Frequency offset in MHz, w.r.t the sequence frequency,
+det = epr.Detection(tp=32,
+                    freq=0, # Frequency offset in GHz, w.r.t the sequence frequency,
                     )
 ```                  
 We now need a time axis for our sequence and to add them to the sequence object.
-When a pulse is copied into the sequence using the `add_pulse` method, parameters can be modified allowing the same pulse can be used multiple times with different timings or amplitudes.
+When a pulse is copied into the sequence using the `addPulse` method, parameters can be modified allowing the same pulse can be used multiple times with different timings or amplitudes.
 ```python
 t = epr.Parameter(name='Interpulse Delay',
                   value=400, # Initial interpulse delay in ns
                   step=8,  # Step size in ns
-                  dim=1024  # Number of points,
-                  unit='ns'  # Unit of the parameter
+                  dim=1024,  # Number of points,
+                  unit='ns',  # Unit of the parameter
                   description='Interpulse delay between the pi/2 and pi pulse'
                   )
 
 # Adding the pulses to the sequence
-seq.add_pulse(p90.copy(t=0))
-seq.add_pulse(p180.copy(t=t))
-seq.add_pulse(det.copy(t=2*t))  
+seq.addPulse(p90.copy(t=0))
+seq.addPulse(p180.copy(t=t))
+seq.addPulse(det.copy(t=2*t))  
 
 # Defining the evolution
 seq.evolution([t])
@@ -81,7 +81,7 @@ HE_Seq = epr.HahnEchoRelaxationSequence(
     shots = 20, # Number of shots per point
     start = 400, # Initial interpulse delay in ns
     step = 8,  # Step size in ns
-    dim = 1024  # Number of points
+    dim = 1024,  # Number of points
     pi2_pulse = p90,  # The 90 degree pulse
     pi_pulse = p180  # The 180 degree pulse
 )

pyepr/classes.py

@@ -45,7 +45,10 @@ def __init__(self,config_file:dict=None,log=None) -> None:
         else:
             self.log = log
         self.resonator = None
-        self.amp_nonlinearity = self.config["Spectrometer"]["Bridge"].get('Amplifier Non-Linearity',None)
+        if self.config != {}:
+            self.amp_nonlinearity = self.config["Spectrometer"]["Bridge"].get('Amplifier Non-Linearity',None)
+        else:
+            self.amp_nonlinearity = None
         pass
 
     def connect(self) -> None:
@@ -505,6 +508,10 @@ def __add__(self, __o:object):
                 raise RuntimeError(
                     "Both parameters axis and the array must have the same shape")
 
+    def __radd__(self, __o:object):
+        return self.__add__(__o) 
+
+
     def __sub__(self, __o:object):
 
         if type(__o) is Parameter:
@@ -573,6 +580,9 @@ def __sub__(self, __o:object):
                 raise RuntimeError(
                     "Both parameters axis and the array must have the same shape")
 
+    def __rsub__(self, __o:object):
+        return self.__sub__(__o)
+
     def __mul__(self, __o:object):
         if type(__o) is Parameter:
             if self.unit != __o.unit:

pyepr/dataset.py

@@ -102,7 +102,8 @@ def create_dataset_from_sequence(data, sequence: Sequence,extra_params={}):
     attr.update({'autoDEER_Version':__version__})
     return xr.DataArray(data, dims=dims, coords=coords,attrs=attr)
 
-def create_dataset_from_axes(data, axes, params: dict = {},axes_labels=None):
+def create_dataset_from_axes(data, axes, params: dict = {},extra_coords:dict = None,
+                             axes_labels=None):
     """
     Create an xarray dataset from a numpy array and a list of axes.
 
@@ -129,7 +130,9 @@ def create_dataset_from_axes(data, axes, params: dict = {},axes_labels=None):
     if not isinstance(axes, list):
         axes = [axes]
     coords = {default_labels.pop(0):a for a in axes}
-    params.update({'autoDEER_Version':__version__})
+    if extra_coords is not None:
+        coords.update(extra_coords)
+    params.update({'PyEPR_Version':__version__})
 
     return xr.DataArray(data, dims=dims, coords=coords, attrs=params)
 
@@ -143,7 +146,7 @@ def create_dataset_from_bruker(filepath):
     labels = []
     for i in range(ndims):
         ax_label = default_labels[i]
-        axis_string = params['DESC'][f'{ax_label}UNI']
+        axis_string = params['DESC'].get(f'{ax_label}UNI',"")
         if "'" in axis_string:
             axis_string = axis_string.replace("'", "")
         if axis_string == 'G':
@@ -162,11 +165,12 @@ def create_dataset_from_bruker(filepath):
 
     coords = {labels[i]:(default_labels[i],a) for i,a in enumerate(axes)}
     attr = {}
-    attr['LO'] = float(params['SPL']['MWFQ'])  / 1e9
-    attr['B'] = float(params['SPL']['B0VL']) * 1e4
-    attr['reptime'] = float(params['DSL']['ftEpr']['ShotRepTime'].replace('us',''))
-    attr['nAvgs'] = int(params['DSL']['recorder']['NbScansAcc'])
-    attr['shots'] = int(params['DSL']['ftEpr']['ShotsPLoop'])
+    if 'DSL' in params:
+        attr['LO'] = float(params['SPL']['MWFQ'])  / 1e9
+        attr['B'] = float(params['SPL']['B0VL']) * 1e4
+        attr['reptime'] = float(params['DSL']['ftEpr']['ShotRepTime'].replace('us',''))
+        attr['nAvgs'] = int(params['DSL']['recorder']['NbScansAcc'])
+        attr['shots'] = int(params['DSL']['ftEpr']['ShotsPLoop'])
     attr.update({'autoDEER_Version':__version__})
 
     return xr.DataArray(data, dims=dims, coords=coords, attrs=attr)
@@ -183,7 +187,7 @@ def save(self, filename, type='netCDF',overwrite=True):
 
         Parameters
         ----------
-        filename : str
+        filename : str, file-like object
             The name of the file to save the dataset
         type : str, optional
             The type of file to save, by default 'netCDF' (including .h5)
@@ -196,7 +200,7 @@ def save(self, filename, type='netCDF',overwrite=True):
             mode = "a"
 
         #if filename doesn't have the extension .h5, add it
-        if not filename.endswith('.h5'):
+        if isinstance(filename,str) and not filename.endswith('.h5'):
             filename = filename + '.h5'
 
         if 'Scan' in self._obj.dims:
@@ -375,3 +379,108 @@ def merge(self,other,ignore_errors=True):
 
 
 
+def downconvert_dataset(dataset, filter_type='boxcar',IF=None,reduce=True,sampling_rate=None,**kwargs):
+    """
+    Downconvert a dataset to baseband using a filter	
+    Parameters
+    ----------
+    dataset : xr.DataArray
+        The dataset to downconvert
+    filter_type : str, optional
+        The type of filter to use, by default 'boxcar'. Other options are 'cheby' and a Pulse object
+    filter_width : float, optional 
+        The width of the filter in MHz or ns depending on filter tyre, by default 20 ns for boxcar and 50 MHz for cheby.
+    IF : float, optional
+        The intermediate frequency to use, by default 0.15
+    reduce : bool, optional
+        If True, the dataset is reduced to a single point, by default True
+    **kwargs : dict
+        Extra arguments to pass to the filter function
+
+    Returns
+    -------
+    xr.DataArray
+        The downconverted dataset
+    """
+
+    if IF is None:
+        if 'IF_freq' in dataset.attrs:
+            IF = dataset.attrs.get('IF_freq') # GHz
+        elif 'if_freq' in dataset.attrs:
+            IF = dataset.attrs.get('if_freq') # GHz
+        elif 'IFfreq' in dataset.attrs:
+            IF = dataset.attrs.get('IFfreq') # GHz
+        else:
+            raise ValueError('IFfreq not found in dataset attributes, please provide IF value')
+
+    if sampling_rate is None:
+        if 'det_rate' in dataset.attrs:
+            sampling_rate = dataset.attrs.get('det_rate') # GHz
+        elif 'sampling_rate' in dataset.attrs:
+            sampling_rate = dataset.attrs.get('sampling_rate') # GHz
+        else:
+            raise ValueError('sampling_rate not found in dataset attributes, please provide sampling_rate value')
+    if filter_type is None:
+        dc_first=True
+    elif filter_type not in ['boxcar','cheby']:
+        filter_type = 'cheby'
+        filter_width = 50
+
+
+    if filter_type == 'boxcar':
+        filter_width = kwargs.get('filter_width',20)
+        funct = lambda data: np.convolve(data,np.ones(filter_width),mode='same')
+        dc_first=True
+    elif isinstance(filter_type, ad_pulses.Pulse): # match filter to a epr Pulse
+        raise NotImplementedError('Match filter to a pulse not implemented yet')
+    elif filter_type == 'cheby':
+        from scipy.signal import cheby1,sosfiltfilt
+        filter_width = kwargs.get('filter_width',50) # MHz
+        filter_width /= 1e3
+
+        Wn = np.array([IF-filter_width,IF+filter_width])
+        Wn[Wn<=0] = 0.001
+        order = 5
+
+        a = cheby1(order,0.5,Wn,'bandpass',analog=False,fs=sampling_rate,output='sos')
+        funct = lambda data: sosfiltfilt(a,data)
+        dc_first=False
+    elif filter_type is None:
+        dc_first=True
+
+    else:
+        raise ValueError('Filter not recognised')
+
+    if dc_first:
+        data_array_dc = dataset*np.exp(-1j*2*np.pi*IF*dataset.tx/sampling_rate)
+        if filter_type is None:
+            return data_array_dc
+        data_array_dc.data = np.apply_along_axis(funct,-1,data_array_dc.data)
+    else:
+        data_array_dc = xr.apply_ufunc(funct,dataset)
+        data_array_dc = data_array_dc*np.exp(-1j*2*np.pi*IF*data_array_dc.tx/sampling_rate)
+
+    # if data_array_dc.ndim == 3:
+    #     max_echo_pos = np.unravel_index(np.argmax(np.abs(data_array_dc.data[0,0,20:-20])),data_array_dc.shape[-1])[0] + 20
+    # else:
+    #     max_echo_pos = np.unravel_index(np.argmax(np.abs(data_array_dc.data[0,20:-20])),data_array_dc.shape[-1])[0] + 20
+    if reduce:
+        if 'max_echo_pos' in kwargs:
+            max_echo_pos = kwargs['max_echo_pos']
+        else:
+            max_echo_pos = np.unravel_index(np.abs(data_array_dc).argmax(),data_array_dc.shape)[-1]
+        return data_array_dc[...,max_echo_pos]
+    else:
+        return data_array_dc
+
+def find_peak(dataset, freq,freq_axis=None,search_range=4):
+    if freq/1e3 < freq_axis.min() or freq/1e3 > freq_axis.max():
+        return "N/A", "N/A"
+
+    if freq_axis is None:
+        freq_axis = dataset.offset_frequency.values
+    n_points = dataset.shape[0]
+    loc = np.argmin(np.abs(freq_axis-freq/1e3))
+    loc = loc-search_range + dataset.values[loc-search_range:loc+search_range].argmax()
+    peak = dataset[loc].values
+    return loc,peak