diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 00000000..709cee15
Binary files /dev/null and b/.DS_Store differ
diff --git a/qualang_tools/.DS_Store b/qualang_tools/.DS_Store
new file mode 100644
index 00000000..80575349
Binary files /dev/null and b/qualang_tools/.DS_Store differ
diff --git a/qualang_tools/config/README_IIR_filters.md b/qualang_tools/config/README_IIR_filters.md
new file mode 100644
index 00000000..f28511e3
--- /dev/null
+++ b/qualang_tools/config/README_IIR_filters.md
@@ -0,0 +1,52 @@
+# Digital filter correction fitting functions
+
+## Overveiw of digital filters
+
+When operating qubits and aiming for high-fidelity single- and two-qubit gates, performance is often limited by the quality of the control pulses generated by the room-temperature controller. Distortions are introduced by the analog setup that connects the controller to the qubits. For example, due to impedance mismatches in the transmission lines. Two main types of filters are typically required:
+
+ - Short time-scale distortions (< 50 ns): These appear as fast, unwanted variations in the control signal. Such deviations can be corrected using FIR (finite impulse response) filters.
+ - Long time-scale distortions (> 50 ns): These manifest as additional RC-like filtering of the signal, leading to long saturation times. Such distortions can be corrected using IIR (infinite impulse response) filters.
+
+This file provides tools for calculating the correction parameters needed to implement digital filters with the OPX1000. In particular, it focuses on generating the IIR filter coefficients that can be directly inserted into the OPX1000 configuration file, making it seamless to implement high-fidelity gates.
+
+More background on digital filters can be found in the [QUA documenation page](https://docs.quantum-machines.co/latest/docs/Guides/output_filter).
+
+## IIR correction filter functions
+
+This repository is optimized for estimating IIR filter corrections and includes tools to fit time-series data with multiple exponential decay components. It is especially useful for analyzing physical systems (e.g., flux or exchange pulse response signals) where the response can be modeled as a sum of exponential decays plus a constant background term.
+
+The workflow assumes that distorted pulse data have already been measured at the qubit level and saved in a ```.h5``` (NetCDF) file.
+
+### Features
+1. Single exponential model implements a simple exponential decay model: f(t) = Ae^(-t/\$\tau$)
+2. Sequential multi-exponential fitting
+ - Automatically estimates a constant (DC) offset from the data tail
+ - Fits the slowest time constants first, subtracts them, and then extracts faster components.
+ - Can constrain decay constants $\tau$ or fit them freely.
+
+### Workflow
+1. Import data in ```.h5``` (NetCDF) file.
+2. Define the ```fitting_start_fractions```. This list defines the number of exponential components (by its size) and provides heuristic, user-chosen starting guesses for their regions.
+These fractions are optimized automatically, but the number of exponents is fixed; in the future, this selection will be automated to avoid user intervention.
+3. Run the optimization ```optimize_start_fractions``` and plot the final result using ```plot_fit```
+5. Check your ```residuals``` to make sure that the fit is reliable
+4. The tuple list of $(A, \tau)$ can now be copied and pasted directly in the configuration file in the ```filter``` -> ```exponential``` section.
+5. Watch the magic unfold and boost your fidelities
+
+### Example
+
+Below is an example of implementing the IIR digital filter using the OPX1000 to correct distortions in a flux-tunable superconducting qubit.
+
+The original uncorrected control pulse was obtained through the [cryoscope experiment](https://pubs.aip.org/aip/apl/article/116/5/054001/38884/Time-domain-characterization-and-correction-of-on), where the LF-FEM channel of the OPX1000 was used to control the qubit flux response. Here the correcsponding [QUA implementatio](https://github.com/qua-platform/qua-libs/blob/main/Quantum-Control-Applications/Superconducting/Single-Flux-Tunable-Transmon/17_cryoscope_1ns.py). Because the flux pulse is distorted by the dilution refrigerator transmission lines, the qubit exhibits a quantum state response directly dependent on the pulse distortion. From this, the actual pulse shape at the qubit can be reconstructed.
+
+As the data are save in the ```.h5``` file, we can now follow the workflow introduced above. In order to succesfully implement the corrections we use the following ```fitting_start_fractions``` = [0.6, 0.3, 0.02] and we obtain the following result
+
+
+
+In this case, only three coefficients were required to reproduce the transfer function of the setup, and these coefficients can be directly used in the configuration file.
+
+Here is a second example, measured on the same setup but from a different qubit. This time we used: ```fitting_start_fractions``` = [0.8, 0.6, 0.3, 0.2, 0.02] and obtained the following result:
+
+
+
+In this second case, five coefficients were required to fully reproduce the transfer function because the high-frequency coaxial cable in the fridge was different.
diff --git a/qualang_tools/config/digital_filters_IIR.py b/qualang_tools/config/digital_filters_IIR.py
new file mode 100644
index 00000000..ab0f4464
--- /dev/null
+++ b/qualang_tools/config/digital_filters_IIR.py
@@ -0,0 +1,293 @@
+# %%
+import numpy as np
+import matplotlib.pyplot as plt
+import netCDF4 as nc
+from typing import List, Tuple
+from scipy.optimize import minimize
+from scipy.optimize import curve_fit
+
+#%%
+def single_exp_decay(t: np.ndarray, amp: float, tau: float) -> np.ndarray:
+ """Single exponential decay without offset
+
+ Args:
+ t (array): Time points
+ amp (float): Amplitude of the decay
+ tau (float): Time constant of the decay
+
+ Returns:
+ array: Exponential decay values
+ """
+ return amp * np.exp(-t/tau)
+
+
+def sequential_exp_fit(
+ t: np.ndarray,
+ y: np.ndarray,
+ start_fractions: List[float],
+ fixed_taus: List[float]=None,
+ a_dc: float=None,
+ verbose: bool=True,
+ ) -> Tuple[List[Tuple[float, float]], float, np.ndarray]:
+ """
+ Fit multiple exponentials sequentially by:
+ 1. First fit a constant term from the tail of the data
+ 2. Fit the longest time constant using the latter part of the data
+ 3. Subtract the fit
+ 4. Repeat for faster components
+
+ Args:
+ t (array): Time points in nanoseconds, representing the time resolution of the pulse.
+ y (array): Amplitude values of the pulse in volts.
+ start_fractions (list): List of fractions (0 to 1) indicating where to start fitting each component. Choice is user defined.
+ fixed_taus (list, optional): Fixed tau values (in nanoseconds) for each exponential component.
+ If provided, only amplitudes are fitted, taus are constrained.
+ Must have same length as start_fractions.
+ a_dc (float, optional): Fixed constant term. If provided, the constant term is not fitted.
+ verbose (bool): Whether to print detailed fitting information
+
+ Returns:
+ tuple: (components, a_dc, residual) where:
+ - components: List of (amplitude, tau) pairs for each fitted component
+ - a_dc: Fitted constant term or the fixed constant term
+ - residual: The difference between the measured data and the fitted curve after subtracting all exponential components.
+ """
+
+ components = [] # List to store (amplitude, tau) pairs
+ t_offset = t - t[0] # Make time start at 0
+
+ # Find the flat region in the tail by looking at local variance
+ window = max(5, len(y) // 20) # Window size by dividing signal into 20 equal pieces or at least 5 points
+ rolling_var = np.array([np.var(y[i:i+window]) for i in range(len(y)-window)])
+ # Find where variance drops below threshold, indicating flat region
+ var_threshold = np.mean(rolling_var) * 0.1 # 10% of mean variance
+ if a_dc is None:
+ try:
+ flat_start = np.where(rolling_var < var_threshold)[0][-1]
+ # Use the flat region to estimate constant term
+ a_dc = np.mean(y[flat_start:])
+ except IndexError:
+ print("No flat region found, using last point of the signal as constant term")
+ a_dc = y[-1]
+
+ if verbose:
+ print(f"\nFitted constant term: {a_dc:.3e}")
+
+ y_residual = y.copy() - a_dc
+
+ for i, start_frac in enumerate(start_fractions):
+ # Calculate start index for this component
+ start_idx = int(len(t) * start_frac)
+ if verbose:
+ print(f"\nFitting component {i+1} using data from t = {t[start_idx]:.1f} ns (fraction: {start_frac:.3f})")
+
+ # Fit current component
+ try:
+ # Prepare fitting parameters based on whether tau is fixed
+ if fixed_taus is not None:
+ # Use fixed tau - only fit amplitude using lambda
+ tau_fixed = fixed_taus[i]
+ p0 = [y_residual[start_idx]] # Only amplitude initial guess
+ if verbose:
+ print(f"Using fixed tau = {tau_fixed:.3f} ns")
+
+ # Perform the fit on the current interval
+ t_fit = t_offset[start_idx:]
+ y_fit = y_residual[start_idx:]
+ popt, _ = curve_fit(lambda t, amp: single_exp_decay(t, amp, tau_fixed), t_fit, y_fit, p0=p0)
+
+ # Store the components
+ amp = popt[0]
+ tau = tau_fixed
+ else:
+ # Fit both amplitude and tau (original behavior)
+ p0 = [
+ y_residual[start_idx], # amplitude
+ t_offset[start_idx] / 3 # tau
+ ]
+
+ # Set bounds for the fit
+ bounds = (
+ [-np.inf, 0.1], # lower bounds: amplitude can be negative, tau must be positive (0.1 ns is arbitrary)
+ [np.inf, np.inf] # upper bounds
+ )
+
+ # Perform the fit on the current interval
+ t_fit = t_offset[start_idx:]
+ y_fit = y_residual[start_idx:]
+ popt, _ = curve_fit(single_exp_decay, t_fit, y_fit, p0=p0, bounds=bounds)
+
+ # Store the components
+ amp, tau = popt
+
+ components.append((amp, tau))
+ if verbose:
+ tau_status = "(fixed)" if fixed_taus is not None else ""
+ print(f"Found component: amplitude = {amp:.3e}, tau = {tau:.3f} ns {tau_status}")
+
+ # Subtract this component from the entire signal
+ y_residual -= amp * np.exp(-t_offset/tau)
+
+ except (RuntimeError, ValueError) as e:
+ if verbose:
+ print(f"Warning: Fitting failed for component {i+1}: {e}")
+ break
+
+ return components, a_dc, y_residual
+
+
+def optimize_start_fractions(t, y, start_fractions, bounds_scale=0.5, fixed_taus=None, a_dc=None, verbose=True):
+ """
+ Optimize the start_fractions by minimizing the RMS between the data and the fitted sum
+ of exponentials using scipy.optimize.minimize.
+
+ Args:
+ t (array): Time points in nanoseconds, representing the time resolution of the pulse.
+ y (array): Amplitude values of the pulse in volts.
+ start_fractions (list): Initial guess for start fractions. Choice is user defined.
+ bounds_scale (float): Scale factor for bounds around start fractions (0.5 means ±50%)
+ fixed_taus (list, optional): Fixed tau values (in nanoseconds) for each exponential component.
+ If provided, only amplitudes are fitted, taus are constrained.
+ Must have same length as start_fractions.
+ a_dc (float, optional): Constant term. If not provided, the constant term is fitted from
+ the tail of the data.
+ verbose (bool): Whether to print detailed fitting information
+
+ Returns:
+ tuple: (success, best_fractions, best_components, best_dc, best_rms)
+ """
+ # Validate fixed_taus parameter
+ if fixed_taus is not None:
+ if len(fixed_taus) != len(start_fractions):
+ raise ValueError("fixed_taus must have the same length as start_fractions")
+ if any(tau <= 0 for tau in fixed_taus):
+ raise ValueError("All fixed_taus values must be positive")
+
+ def objective(x):
+ """
+ Objective function to minimize: RMS between the data and the fitted sum of
+ exponentials.
+ """
+ # Ensure fractions are ordered in descending order
+ if not np.all(np.diff(x) < 0):
+ return 1e6 # Return large value if constraint is violated
+
+ components, _, residual = sequential_exp_fit(t, y, x, fixed_taus=fixed_taus, a_dc=a_dc, verbose=verbose)
+ if len(components) == len(start_fractions):
+ current_rms = np.sqrt(np.mean(residual**2))
+ else:
+ current_rms = 1e6 # Return large value if fitting fails
+
+ return current_rms
+
+
+ # Define bounds for optimization
+ bounds = []
+ for start in start_fractions:
+ min_val = start * (1 - bounds_scale)
+ max_val = start * (1 + bounds_scale)
+ bounds.append((min_val, max_val))
+
+ print("\nOptimizing start_fractions using scipy.optimize.minimize...")
+ print(f"Initial values: {[f'{f:.5f}' for f in start_fractions]}")
+ print(f"Bounds: ±{bounds_scale*100}% around initial values")
+
+ # Run optimization
+ result = minimize(
+ objective,
+ x0=start_fractions,
+ bounds=bounds,
+ method='Nelder-Mead', # This method works well for non-smooth functions
+ options={'disp': True, 'maxiter': 200}
+ )
+
+ # Get final results
+ if result.success:
+ best_fractions = result.x
+ components, a_dc, best_residual = sequential_exp_fit(t, y, best_fractions, fixed_taus=fixed_taus, a_dc=a_dc, verbose=False)
+ best_rms = np.sqrt(np.mean(best_residual**2))
+ print("\nOptimization successful!")
+ print(f"Initial fractions: {[f'{f:.5f}' for f in start_fractions]}")
+ print(f"Optimized fractions: {[f'{f:.5f}' for f in best_fractions]}")
+ if fixed_taus is not None:
+ print(f"Fixed taus: {[f'{tau:.3f} ns' for tau in fixed_taus]}")
+ print(f"Final RMS: {best_rms:.3e}")
+ print(f"Number of iterations: {result.nit}")
+ else:
+ print("\nOptimization failed. Using initial values.")
+ best_fractions = start_fractions
+ components, a_dc, best_residual = sequential_exp_fit(t, y, best_fractions, fixed_taus=fixed_taus, a_dc=a_dc, verbose=False)
+ best_rms = np.sqrt(np.mean(best_residual**2))
+
+ components = [(amp * np.exp(t[0] / tau), tau) for amp, tau in components]
+ print(components)
+ return result.success, best_fractions, components, a_dc, best_rms
+
+def plot_fit(t_data: np.ndarray, y_data: np.ndarray, components: List[Tuple[float, float]], a_dc: float):
+ """Plot exponential fit results with both linear and log scales.
+
+ Args:
+ t_data (np.ndarray): Time points in nanoseconds
+ y_data (np.ndarray): Measured flux response data
+ components (List[Tuple[float, float]]): List of (amplitude, tau) pairs for each fitted component
+ a_dc (float): Constant term
+
+ Returns:
+ tuple: (fig, axs) where:
+ - fig: Figure object
+ - axs: List of axes objects
+ """
+
+ fit_text = f'a_dc = {a_dc:.3f}\n'
+ y_fit = np.ones_like(t_data, dtype=float) * a_dc
+ for i, (amp, tau) in enumerate(components):
+ y_fit += amp * np.exp(-t_data/tau)
+ fit_text += f'a{i+1} = {amp / a_dc:.3f}, τ{i+1} = {tau:.0f}ns\n'
+
+ fig, axs = plt.subplots(1, 2, figsize=(12, 5))
+
+ # First subplot - linear scale
+ axs[0].plot(t_data, y_data, label='Data')
+ axs[0].plot(t_data, y_fit, label='Fit')
+ axs[0].text(0.98, 0.5, fit_text, transform=axs[0].transAxes, fontsize=10,
+ horizontalalignment='right', verticalalignment='center')
+ axs[0].set_xlabel('Time (ns)')
+ axs[0].set_ylabel('Flux Response')
+ axs[0].legend()
+ axs[0].grid(True)
+ axs[0].ticklabel_format(axis='x', style='sci', scilimits=(0,0))
+
+ # Second subplot - log scale
+ axs[1].plot(t_data, y_data, label='Data')
+ axs[1].plot(t_data, y_fit, label='Fit')
+ axs[1].text(0.98, 0.5, fit_text, transform=axs[1].transAxes, fontsize=10,
+ horizontalalignment='right', verticalalignment='center')
+ axs[1].set_xlabel('Time (ns)')
+ axs[1].set_ylabel('Flux Response')
+ axs[1].set_xscale('log')
+ axs[1].legend(loc='best')
+ axs[1].grid(True)
+
+ fig.tight_layout()
+
+ return fig, axs
+
+# %%
+# load data
+with nc.Dataset('/Users/fabioansaloni/Downloads/ds2.h5', 'r') as f:
+ t_data = np.array(f.variables['time'][:])
+ y_data = np.squeeze(np.array(f.variables['flux_response'][:]))
+
+# fit data
+fitting_start_fractions = [0.8, 0.6, 0.3, 0.2, 0.02]
+# fitting_start_fractions = [0.6, 0.3, 0.02]
+success, best_fractions, components, a_dc, best_rms = optimize_start_fractions(t_data, y_data,fitting_start_fractions)
+
+# plot data and fit
+fig, axs = plot_fit(t_data, y_data, components, a_dc)
+
+# parameters to update in config
+A_list = [component[0] / a_dc for component in components]
+tau_list = [component[1] for component in components]
+exponential_filter = list(zip(A_list, tau_list))
+# %%
diff --git a/qualang_tools/config/image.png b/qualang_tools/config/image.png
new file mode 100644
index 00000000..d81449f4
Binary files /dev/null and b/qualang_tools/config/image.png differ
diff --git a/qualang_tools/config/image1.png b/qualang_tools/config/image1.png
new file mode 100644
index 00000000..16c48b46
Binary files /dev/null and b/qualang_tools/config/image1.png differ
diff --git a/qualang_tools/digital_filters/FIR_first.png b/qualang_tools/digital_filters/FIR_first.png
new file mode 100644
index 00000000..c9154671
Binary files /dev/null and b/qualang_tools/digital_filters/FIR_first.png differ
diff --git a/qualang_tools/digital_filters/FIR_second.png b/qualang_tools/digital_filters/FIR_second.png
new file mode 100644
index 00000000..1d069672
Binary files /dev/null and b/qualang_tools/digital_filters/FIR_second.png differ
diff --git a/qualang_tools/digital_filters/README_FIR_filters.md b/qualang_tools/digital_filters/README_FIR_filters.md
new file mode 100644
index 00000000..bb385557
--- /dev/null
+++ b/qualang_tools/digital_filters/README_FIR_filters.md
@@ -0,0 +1,118 @@
+
+# Digital filter correction fitting functions
+
+## Overview of digital filters
+
+When operating qubits and aiming for high-fidelity single- and two-qubit gates, performance is often limited by the quality of the control pulses generated by the room-temperature controller. Distortions are introduced by the analog setup that connects the controller to the qubits. For example, due to impedance mismatches in the transmission lines. Two main types of filters are typically required:
+
+ - Short timescale distortions ($\lesssim$ 30 ns): These appear as fast, unwanted variations in the control signal. Such deviations can be corrected using FIR (finite impulse response) filters.
+ - Long timescale distortions ($\gtrsim$ 30 ns): These manifest as additional RC-like filtering of the signal, leading to long saturation times. Such distortions can be corrected using IIR (infinite impulse response) filters.
+
+This file provides tools for calculating the correction parameters needed to implement digital filters with the OPX1000. In particular, it focuses on generating the IIR filter coefficients that can be directly inserted into the OPX1000 configuration file, making it seamless to implement high-fidelity gates.
+
+More background on digital filters can be found in the [QUA documentation page](https://docs.quantum-machines.co/latest/docs/Guides/output_filter).
+
+## FIR correction filter functions
+
+This repository is optimized for estimating FIR filter corrections and it is optimized to work once IIR corrections have already been extracted. The following FIR implemetnation is based on the procedure described in this [paper](https://arxiv.org/abs/2503.04610).
+
+The FIR filters works at 2 GSa/s, and there are 48 taps availbale per core, therefoe allowing corrections for the first 24 ns. \
+First, we extract the forward FIR filter that, when convolved with an ideal step, best reproduces the measured response; this is done by solving a regularized least-squares problem that balances reconstruction accuracy with smooth, physically plausible coefficients. A Toeplitz convolution matrix is used to express the convolution. Noise overfitting is prevetned by implementing a Tikhonov regulizer while a second regulizer is used to physically model that the filter deos not extend indefinitely.\
+The best combination of FIR length and regularization parameters is selected by minimizing the reconstruction error across candidate models. \
+After the forward FIR is found, a second optimization computes an inverse FIR filter that, when applied before the hardware response, compensates distortions by approximating a delta function. A Gaussian target and a smoothness regularizer make this inversion stable.
+
+The workflow assumes that distorted pulse data have already been measured at the qubit level and corrected for IIR distorsions. After that data should be saved in a ```.h5``` (NetCDF) file.
+
+### Features
+1. Forward FIR system identification
+ - Fits a finite-impulse-response (FIR) filter that transforms an ideal step into the measured step response using regularized least-squares.
+ - Uses multiple regularization coefficients to suppress noise amplifications.
+ - Automatically optimizes FIR lengths and regularizations parameters.
+2. Inverse FIR filters extraction - the predistorsion filters
+ - Computes a second FIR filter whose convolution with the forward filter approximates a delta function.
+ - Uses a Gaussian target as surrogate for the ideal delta and a smoothness regularizer to avoid unstable solutions from the inversion procedure
+
+### Workflow
+1. Import data in ```.h5``` (NetCDF) file. Remeber to correct the data for IIR distorsion before first.
+2. Set the ```FIR_filters_coefficitents```: Tikhonow regularization parameters ```lam1_values```, exponential decay regularization parameters ```lam2_values``` as well as the number of FIR coefficients ```L```. These coefficients try to model the setup transfer function. Since this can be complicated, the algorithm evaluates many possible FIR filters and chooses the one that best explains the measured data.
+3. Run the optimization ```analyze_and_plot_fir_fit``` (this is automatically called in the main function). For each parameter combination, the function reconstructs the step response and chooses the combination that gives the best fit. The function plots the measured vs reconstructed step response and prints the residuals as well as the stability of the FIR coefficients. If the fit is smooth an residual are small, continue to next step. If results are inconsisitend, manually optimize the FIR filter coefficients.
+4. Now set the ```inverse_FIR_filter_coefficients```: the ```method``` to use for inversion of the previusly calculated transfer function (default is ```optimization```), the ```sigma_ns``` which defines the width of the Gaussian which approximates the delta function resulting from the inversion (typically 0.5ns or 0.25ns is used) and eventually ```lam_smooth``` is the regularization parameter used to limit rapid changes in the inverse FIR coefficients. Smaller values of ```sigma_ns``` and ```lam_smooth``` produce more aggressive correction at the cost of noise amplification and possible instabilities. Values in the range ```sigma_ns``` = 0.25-0.75ns and ```lam_smooth``` = 0.01 - 0.1 should suffice but the user can explore other possibilities based on the data at hand.
+5. Run the optimization ```analyze_and_plot_inverse_fir```, which for each combination of ```inverse_FIR_filter_coefficients``` it solves the inverse step response. The function plots the measured vs reconstructed step response and displays the residuals along with the stability of the FIR coefficients. If the reconstruction is smooth, residuals are small, and the predistortion looks physically reasonable, the FIR fit is considered valid and you can proceed. If the results are inconsistent, manually adjust the ```inverse_FIR_filter_coefficients``` and rerun the analysis.
+6. The list of ```fir_filter``` can now be copied and pasted directly in the configuration file in the ```filter``` -> ```feedforward``` section.
+7. Watch the magic unfold and boost your gate fidelity even more!
+
+### Example
+
+Below is an example of implementing the FIR digital filter using the OPX1000 to correct distortions in a flux-tunable superconducting qubit.
+
+The original uncorrected control pulse was obtained through the [cryoscope experiment](https://pubs.aip.org/aip/apl/article/116/5/054001/38884/Time-domain-characterization-and-correction-of-on), where the LF-FEM channel of the OPX1000 was used to control the qubit flux response. Here the corresponding [QUA implementation](https://github.com/qua-platform/qua-libs/blob/main/Quantum-Control-Applications/Superconducting/Single-Flux-Tunable-Transmon/17_cryoscope_1ns.py). Because the flux pulse is distorted by the dilution refrigerator transmission lines, the qubit exhibits a quantum state response directly dependent on the pulse distortion. From this, the actual pulse shape at the qubit can be reconstructed.
+
+As the data are saved in the ```.h5``` file, we can now follow the workflow introduced above:
+````python
+import numpy as np
+import netCDF4 as nc
+
+# load data
+with nc.Dataset('path_to_dataset.h5', 'r') as f:
+ t_data = np.array(f.variables['time'][:])
+ y_data = np.squeeze(np.array(f.variables['flux_response'][:]))
+
+# normalize the response and resample the response to 2 GS/s
+normalized_response_raw = y_data / y_data[-10:].mean()
+original_Ts = t_data[1] - t_data[0]
+normalized_response_2gsps = resample_to_target_rate(normalized_response_raw, original_Ts, 0.5)
+time_2gsps = np.arange(len(normalized_response_2gsps)) * 0.5 + 1
+
+#define the coefficients
+FIR_filter_coefficient = [[16, 20, 24, 28, 32, 40, 48],[1e-5, 1e-4, 1e-3, 1e-2, 1e-1],[1e-5, 1e-4, 1e-3, 1e-2, 1e-1]]
+inverse_FIR_filter_coefficient =[0.75,5e-2,'optimization' ]
+
+# fit the data and plot the forward and inverse FIR filters
+h_fir, inv_fir, fig_fir_fit, fig_inv_fir_fit = analyze_and_plot_inverse_fir(
+ response=normalized_response_2gsps,
+ time=time_2gsps
+ )
+
+# calculate the predistorted and corrected response
+ideal_response = np.ones(len(y_data))
+predistorted_response = lfilter(inv_fir, 1, ideal_response)
+corrected_response = lfilter(h_fir, 1, predistorted_response)
+
+# feedforward filters to be used in config file
+fir_filter = inv_fir.tolist()
+
+# Assuming that we want to update the filters on analog output port 1 of LF-fem 1
+output_port = config["controllers"]["con1"]["fems"][1]["analog_outputs"][1]
+output_port["filter"] = {
+ "exponential_dc_gain": a_dc,
+ "exponential": exponential_filter,
+ "feedforward": fir_filter
+}
+
+
+````
+On a given qubit, after the IIR filter function has been applied, we would obtain a step response at the qubit that will look like the following one
+
+
+
+By selecting the ```FIR_filters_coefficitents``` = [[16, 20, 24, 28, 32, 40, 48],[1e-5, 1e-4, 1e-3, 1e-2, 1e-1],[1e-5, 1e-4, 1e-3, 1e-2, 1e-1]] (where the coefficients are lists so that the algorithm can iterate over multiple configurations), the result of the calculated hardware setup transfer function will look like the follwoing plot:
+
+
+
+The upper panels show the reconstructed and measured normalized step response of Fig. 1 (left panel) and the correpsonding difference (right panel). The bottom right panel shows the mean and standard deviation of the rms reconstruction error for each number of FIR coefficients, with the lowest error in red horizontal dashed line. The best FIR coefficients are shown on the bottom left panel. This fit is succesful with BEST PARAMETERS =[L = 48, lam1 = 1e-03, lam2 = 1e-04] and a reconstruction error of 1.6763e-03 (0.17%).
+
+We can now move forward and define the ```inverse_FIR_filter_coefficients``` = [0.75,5e-2,'optimization' ] and proceed to extract the inverse of the hardware setup response which will be used in the FIR filters coefficients in the OPX100 configuration file. For these coefficients, we obtain the following results:
+
+
+
+The upper left panel shows the measured and reconstructed signal based on the forward FIR extraction (similar to the upper left panel of Fig. 3), and upper right panel shows the extracted forward and inverse FIR coefficients. The middle left panel shows the predistorted signal and the expected corrected signal, and the middle right panel shows also the measured signal and the backward corrected signal (calculated from the measured signal and the inverse FIR filter). The lower panels show the difference between the corrected signal and the ideal step function (left) and the convolution of the forward and inverse FIR filters, that should be close to the target Gaussian approximating a delta function. The obtained correction error in the is 5.713e-02 and therefore acceptable.
+
+The FIR filter coeffieicent are now saved in the configuration file. By now performing a new cryoscope experimetn (where the predistorted pulse is played), we obtain the following mesaured step response at the qubit (in blue)
+
+
+
+One can perform the FIR filter estimation again, whihc would convolute the existing filters and the new extracted ones. Such iterations can be repeated until no more improvement is observed. We eventually achieved the following result
+
+
+
+Showing an effective rise time of <1 ns (as the first point of 1ns is set at 0.97) as well as a clear saturation within 0.1% of the signal within 5 ns. In both figures, the black solid lines represent 0.01% of the target step impulse (normalized to 1).
diff --git a/qualang_tools/digital_filters/fir_filters_fit.py b/qualang_tools/digital_filters/fir_filters_fit.py
new file mode 100644
index 00000000..6dec48cd
--- /dev/null
+++ b/qualang_tools/digital_filters/fir_filters_fit.py
@@ -0,0 +1,746 @@
+import matplotlib.pylab as plt
+import numpy as np
+from scipy import linalg
+from scipy.interpolate import interp1d
+from scipy.signal import lfilter
+from typing import List, Tuple, Literal
+from typing import Optional
+from numpy.typing import ArrayLike
+import netCDF4 as nc
+
+#%%
+def conv_causal(v: ArrayLike, h: ArrayLike, N: Optional[int] = None) -> np.ndarray:
+ """
+ Perform a causal (one-sided) convolution of input signal v with filter h.
+
+ Parameters
+ ----------
+ v : array-like
+ Input sequence (e.g., signal to be filtered)
+ h : array-like
+ Impulse response (filter coefficients)
+ N : int or None, optional
+ Number of output points to return. If None, returns the convolution up to the length of v.
+
+ Returns
+ -------
+ y : ndarray
+ The result of the causal convolution of v with h, truncated to N or len(v) samples.
+ """
+ v = np.asarray(v, dtype=float)
+ h = np.asarray(h, dtype=float)
+ y = np.convolve(v, h, mode='full')
+
+ return y[: (len(v) if N is None else N)]
+
+
+def build_toeplitz_matrix(v, L):
+ """
+ Construct a Toeplitz matrix from input sequence v for FIR system identification.
+
+ Parameters
+ ----------
+ v : array-like
+ Input sequence (stimulus to system).
+ L : int
+ Number of filter coefficients (FIR length).
+
+ Returns
+ -------
+ V : ndarray, shape (N, L)
+ Toeplitz matrix such that each row forms a lagged vector of v,
+ suitable for linear convolution: phi ≈ V @ h.
+ """
+ v = np.asarray(v, float)
+ N = len(v)
+
+ return linalg.toeplitz(c=v, r=np.concatenate([[v[0]], np.zeros(L - 1)]))
+
+
+def resample_to_target_rate(
+ data: ArrayLike,
+ original_Ts: float,
+ target_Ts: float,
+ kind: str = 'cubic'
+ ) -> np.ndarray:
+ """
+ Resample time-domain data to a specified target sampling rate.
+
+ Parameters
+ ----------
+ data : array-like
+ The original time-domain data
+ original_Ts : float
+ The original sampling time, in nanoseconds
+ target_Ts : float
+ The target sampling time, in nanoseconds
+ kind : str
+ Interpolation kind ('linear', 'cubic', etc.)
+
+ Returns
+ -------
+ resampled : numpy.ndarray
+ The data interpolated onto the target sampling rate grid
+ """
+
+ data = np.asarray(data)
+ N = len(data)
+ t_original = np.arange(N) * original_Ts
+ max_time = t_original[-1]
+
+ num_samples = int(np.floor(max_time / target_Ts)) + 1
+ t_target = np.arange(num_samples) * target_Ts
+ t_target = t_target[t_target <= max_time]
+
+ interp_fun = interp1d(t_original, data, kind=kind, fill_value="extrapolate", bounds_error=False)
+ return interp_fun(t_target)
+
+FIR_filter_coefficient =[]
+
+def fit_fir(
+ phi: ArrayLike,
+ v: ArrayLike,
+ L: int = 0.5,
+ Ts: float = 0.5,
+ lam1: float = 1e-2,
+ lam2: float = 1e-2,
+ tail_ns: Optional[float] = None
+ ) -> np.ndarray:
+ """
+ Fit a finite impulse response (FIR) filter to data. Solves for FIR filter coefficients `h`
+ with length `L` such that `phi ≈ V @ h`, where `V` is the Toeplitz matrix constructed from
+ input `v`.
+
+ Parameters
+ ----------
+ phi : array_like
+ response signal (1D array).
+ v : array_like
+ Input signal to be filtered (1D array).
+ L : int
+ Number of FIR filter taps (length of FIR filter).
+ Ts : float, optional
+ Sampling time (in nanoseconds). Default is 0.5.
+ lam1 : float, optional
+ Regularization parameter for the identity matrix. Default is 1e-2.
+ lam2 : float, optional
+ Regularization parameter for exponential tail. Default is 1e-2.
+ tail_ns : float, optional
+ Time constant (in nanoseconds) for tail regularization. If None, computed as (L*Ts)/3.0.
+
+ Returns
+ -------
+ h : ndarray
+ Fitted FIR filter coefficients (1D array of length L).
+ """
+ phi = np.asarray(phi, float)
+ v = np.asarray(v, float)
+ V = build_toeplitz_matrix(v, L)
+
+ if tail_ns is None:
+ tail_ns = (L*Ts)/3.0
+ idx = np.arange(L)
+ x = np.exp(idx * Ts / tail_ns)
+
+ A = V.T @ V + lam1*np.eye(L) + lam2*np.diag(x)
+ b = V.T @ phi
+ h = linalg.solve(A, b, assume_a='pos')
+
+ return h
+
+
+def optimize_fir_parameters(
+ response: ArrayLike,
+ Ts: float = 0.5,
+ L_values: List[int] = [16, 20, 24, 28, 32, 40, 48],
+ lam1_values: List[float] = [1e-5, 1e-4, 1e-3, 1e-2, 1e-1],
+ lam2_values: List[float] = [1e-5, 1e-4, 1e-3, 1e-2, 1e-1]
+ ) -> Tuple[List[dict], float, dict, np.ndarray, np.ndarray]:
+ """
+ Optimize FIR extraction parameters (L, lam1, lam2) to minimize the reconstruction error
+ between the measured response and the reconstructed signal.
+
+ Parameters
+ ----------
+ response: array_like
+ response signal (1D array).
+ Ts: float, optional
+ Sampling time (in nanoseconds). Default is 0.5.
+ L_values: List[int], optional
+ List of FIR filter lengths to search over. Default is [16, 20, 24, 28, 32, 40, 48].
+ lam1_values: List[float], optional
+ List of regularization parameters for the identity matrix to search over. Default is [1e-5, 1e-4, 1e-3, 1e-2, 1e-1].
+ lam2_values: List[float], optional
+ List of regularization parameters for the exponential tail to search over. Default is [1e-5, 1e-4, 1e-3, 1e-2, 1e-1].
+
+ Returns
+ -------
+ best_error: minimum NRMS error found
+ best_params: dictionary of optimal parameters (L, lam1, lam2, error, h, reconstructed)
+ best_h: best FIR filter coefficients found
+ best_reconstructed: reconstructed signal using the best FIR filter
+ results: list of dictionaries of all results for each parameter combination
+ """
+ print("="*70)
+ print("OPTIMIZING FIR EXTRACTION PARAMETERS")
+ print("="*70)
+ print("Searching over L, lam1, lam2 to minimize reconstruction error")
+ print("Error = ||response - reconstructed|| / ||response||\n")
+ print(f"Signal length: {len(response)} samples")
+ print(f"Time span: {(len(response) * Ts):.1f} ns")
+ print(f"Sampling time: {Ts:.2f} ns")
+ print(f"Sampling rate: {1/Ts:.2f} GS/s\n")
+
+ total_combinations = len(L_values) * len(lam1_values) * len(lam2_values)
+ print(f"Total parameter combinations to test: {total_combinations}")
+ print(f"L: {L_values}")
+ print(f"lam1: {lam1_values}")
+ print(f"lam2: {lam2_values}")
+ print("\nStarting optimization...\n")
+
+ results = []
+ best_error = float('inf')
+ best_params = None
+ best_h = None
+ best_reconstructed = None
+
+ # Ideal signal assumed as a step function
+ ideal_response = np.ones(len(response))
+
+ iteration = 0
+ for L in L_values:
+ for lam1 in lam1_values:
+ for lam2 in lam2_values:
+ iteration += 1
+ try:
+ # Extract forward FIR filter
+ h = fit_fir(response, ideal_response, L=L, Ts=Ts, lam1=lam1, lam2=lam2)
+ h /= np.sum(h) # Normalize
+
+ # Reconstruct signal using Toeplitz matrix
+ V = build_toeplitz_matrix(ideal_response, L)
+ reconstructed = V @ h
+
+ # # Truncate or pad reconstructed to match length
+ # if len(reconstructed) > len(response):
+ # reconstructed = reconstructed[:len(response)]
+ # elif len(reconstructed) < len(response):
+ # reconstructed = np.pad(reconstructed, (0, len(response) - len(reconstructed)),
+ # mode='edge')
+
+ # Compute reconstruction error (NRMS)
+ reconstruction_error = np.linalg.norm(response - reconstructed) / np.linalg.norm(response)
+
+ # Store results
+ result = {
+ 'L': L,
+ 'lam1': lam1,
+ 'lam2': lam2,
+ 'error': reconstruction_error,
+ 'h': h.copy(),
+ 'reconstructed': reconstructed.copy()
+ }
+ results.append(result)
+
+ # Track best
+ if reconstruction_error < best_error:
+ best_error = reconstruction_error
+ best_params = result
+ best_h = h.copy()
+ best_reconstructed = reconstructed.copy()
+
+ except Exception as e:
+ print(f" Warning: Failed for L={L}, lam1={lam1:.0e}, lam2={lam2:.0e}: {e}")
+ continue
+
+ print(f"\n{'='*70}")
+ print(f"OPTIMIZATION COMPLETE")
+ print(f"{'='*70}")
+ print(f"Total combinations tested: {len(results)}")
+
+ return results, best_error, best_params, best_h, best_reconstructed
+
+
+def analyze_and_plot_fir_fit(
+ response: ArrayLike,
+ time: ArrayLike,
+ Ts: float = 0.5,
+ L_values: List[int] = [16, 20, 24, 28, 32, 40, 48],
+ lam1_values: List[float] = [1e-5, 1e-4, 1e-3, 1e-2, 1e-1],
+ lam2_values: List[float] = [1e-5, 1e-4, 1e-3, 1e-2, 1e-1],
+ verbose: bool = True
+ ) -> Tuple[np.ndarray, dict, np.ndarray, float]:
+ """
+ For a given response signal, analyze the FIR filter fit that best reconstructs the response
+ from the ideal step response and plot the results.
+
+ Parameters
+ ----------
+ response : ArrayLike
+ The measured (distorted) response signal to analyze, typically normalized amplitude data.
+ time : ArrayLike
+ The corresponding time array for the response signal (in ns).
+ Ts : float, optional
+ The sampling interval (in ns). Default is 0.5.
+ L_values : List[int], optional
+ The list of FIR filter lengths to search over. Default is [16, 20, 24, 28, 32, 40, 48].
+ lam1_values : List[float], optional
+ The list of regularization parameters for the identity matrix to search over. Default is [1e-5, 1e-4, 1e-3, 1e-2, 1e-1].
+ lam2_values : List[float], optional
+ The list of regularization parameters for the exponential tail to search over. Default is [1e-5, 1e-4, 1e-3, 1e-2, 1e-1].
+ verbose: bool, optional
+ Whether to print verbose output. Default is True.
+
+ Returns
+ -------
+ best_h : np.ndarray
+ The coefficients of the best forward FIR filter found.
+ best_params : dict
+ Dictionary of the optimal FIR filter hyperparameters.
+ best_reconstructed : np.ndarray
+ The best reconstructed (filtered) signal from the optimization.
+ best_error : float
+ Normalized root mean square error (NRMS) for the optimal filter parameters.
+ fig : matplotlib.figure.Figure
+ The figure object containing the signal reconstruction and FIR coefficients plots.
+ """
+ results, best_error, best_params, best_h, best_reconstructed = optimize_fir_parameters(
+ response, Ts=Ts, L_values=L_values, lam1_values=lam1_values, lam2_values=lam2_values
+ )
+ if best_params is not None and verbose:
+ print(f"\nBEST PARAMETERS:")
+ print(f" L = {best_params['L']}")
+ print(f" lam1 = {best_params['lam1']:.0e}")
+ print(f" lam2 = {best_params['lam2']:.0e}")
+ print(f"\nPERFORMANCE:")
+ print(f" Reconstruction error: {best_error:.4e} ({best_error*100:.2f}%)")
+ print(f" Max |h|: {np.max(np.abs(best_h)):.4e}")
+ print(f" h.sum(): {np.sum(best_h):.6f}")
+
+ # Visualize best result
+ fig, axes = plt.subplots(2, 2, figsize=(14, 10))
+
+ # Plot 1: Signal comparison
+ ax = axes[0, 0]
+ ax.plot(time, response, 'r-', label='Measured (distorted)', linewidth=2, alpha=0.7)
+ ax.plot(time, best_reconstructed, 'b--',
+ label=f'Reconstructed (error={best_error:.4e})', linewidth=2)
+ ax.set_xlabel('Time (ns)')
+ ax.set_ylabel('Amplitude')
+ ax.set_title('Best Reconstruction Result')
+ ax.legend()
+ ax.grid(True, alpha=0.3)
+
+ # Plot 2: Residual
+ ax = axes[0, 1]
+ residual = response - best_reconstructed
+ ax.plot(time, residual, 'm-', linewidth=1.5)
+ ax.axhline(y=0, color='k', linestyle='--', alpha=0.3)
+ ax.fill_between(time, -np.std(residual), np.std(residual), alpha=0.2, color='gray')
+ ax.set_xlabel('Time (ns)')
+ ax.set_ylabel('Residual')
+ ax.set_title(f'Reconstruction Residual (σ={np.std(residual):.4e})')
+ ax.grid(True, alpha=0.3)
+
+ # Plot 3: Best FIR filter
+ ax = axes[1, 0]
+ ax.plot(best_h, 'b-o', markersize=4, linewidth=2)
+ ax.set_xlabel('Tap Index')
+ ax.set_ylabel('Coefficient Value')
+ ax.set_title(f'Best Forward FIR (L={best_params["L"]})')
+ ax.grid(True, alpha=0.3)
+
+ # Plot 4: Error distribution across parameters
+ ax = axes[1, 1]
+ errors_by_L = {}
+ for r in results:
+ if r['L'] not in errors_by_L:
+ errors_by_L[r['L']] = []
+ errors_by_L[r['L']].append(r['error'])
+
+ L_sorted = sorted(errors_by_L.keys())
+ error_means = [np.mean(errors_by_L[L]) for L in L_sorted]
+ error_stds = [np.std(errors_by_L[L]) for L in L_sorted]
+
+ ax.errorbar(L_sorted, error_means, yerr=error_stds, fmt='o-', capsize=5, linewidth=2, markersize=6)
+ ax.axhline(y=best_error, color='r', linestyle='--', alpha=0.5, label=f'Best: {best_error:.4e}')
+ ax.set_xlabel('Filter Length L')
+ ax.set_ylabel('Reconstruction Error')
+ ax.set_title('Error vs Filter Length')
+ ax.legend()
+ ax.grid(True, alpha=0.3)
+
+ plt.tight_layout()
+ plt.show()
+
+ # Show top 10 results
+ print(f"\n{'='*70}")
+ print("TOP 10 PARAMETER COMBINATIONS:")
+ print(f"{'='*70}")
+ sorted_results = sorted(results, key=lambda x: x['error'])[:10]
+ print(f"{'Rank':<6} {'L':<4} {'lam1':<8} {'lam2':<8} {'Error':<12}")
+ print("-"*70)
+ for i, r in enumerate(sorted_results, 1):
+ print(f"{i:<6} {r['L']:<4} {r['lam1']:<8.0e} {r['lam2']:<8.0e} {r['error']:<12.4e}")
+ print("="*70)
+
+ elif verbose:
+ print("ERROR: No valid parameter combinations found!")
+ print("Try adjusting parameter ranges or check your data.")
+
+ return best_h, best_params, best_reconstructed, best_error, fig
+
+
+def invert_fir(
+ h: ArrayLike,
+ Ts: float = 0.5,
+ M: Optional[int] = None,
+ method: Literal['optimization', 'analytical'] = 'optimization',
+ sigma_ns: float = 0.75,
+ lam_smooth: float = 5e-2 ,
+ normalize_dc_gain: bool = False
+ ) -> np.ndarray:
+ """
+ Invert a causal FIR filter by finding an inverse causal FIR (possibly smoothed).
+
+ Solves:
+ min_{h_inv} || d - (h * h_inv) ||^2 + lam_smooth * || Δ h_inv ||^2
+
+ where:
+ - h is the original causal FIR filter coefficients (array-like, length L)
+ - h_inv is the inverse FIR filter coefficients (array-like, length M)
+ - d is a desired target response (by default, a normalized Gaussian to approximate a causal impulse)
+ - * denotes convolution (implemented via Toeplitz matrix multiplication)
+ - Δ is the first difference operator (to penalize roughness in h_inv)
+ - lam_smooth is the regularization/smoothing parameter
+
+ Parameters
+ ----------
+ h : array_like
+ Original FIR filter coefficients (causal, length L).
+ Ts : float, optional
+ Sample spacing (default is 0.5).
+ M : int, optional
+ Length of the inverse FIR filter coefficients to solve for (default is len(h)).
+ method : str, optional
+ Method to use for inversion ('optimization' or 'analytical').
+ sigma_ns : float, optional
+ Standard deviation of target Gaussian for delta approximation (default is 1.0).
+ lam_smooth : float, optional
+ Regularization parameter for first-difference smoothing of h_inv (default is 5e-2).
+ normalize_dc_gain: bool, optional
+ Whether to normalize the DC gain of the inverse FIR filter to 1. Default is False.
+
+ Returns
+ -------
+ h_inv : ndarray
+ FIR inverse coefficients (length M), optionally DC gain-normalized.
+
+ Notes
+ -----
+ - This routine finds a stable, causal FIR approximate inverse of a given FIR.
+ - Regularization improves stability for nearly non-minimum-phase FIRs.
+ - The returned inverse may not be exact due to smoothing and length limits.
+ """
+ h = np.asarray(h, float)
+ L = len(h)
+ if M is None:
+ M = L
+ if method == 'optimization':
+ # target 'delta'
+ t = np.arange(M)*Ts
+ d = np.exp(-0.5*(t/(sigma_ns))**2)
+ d /= d.sum() # unit DC gain
+
+ # Build Toeplitz conv matrix H for causal conv(h, h_inv) truncated to M
+ H = build_toeplitz_matrix(h, M)[:M, :]
+
+ # First-difference (Sobolev) smoothing
+ D = np.eye(M, k=0) - np.eye(M, k=1)
+ D = D[:-1, :] # (M-1) x M
+
+ A = H.T @ H + lam_smooth * (D.T @ D)
+ b = H.T @ d
+ h_inv = linalg.solve(A, b, assume_a='pos')
+
+ # Optional: normalize composite DC gain to 1
+ if normalize_dc_gain:
+ gain = (h.sum() * h_inv.sum())
+ if gain != 0:
+ h_inv /= gain
+
+ elif method == 'analytical':
+ h_inv = np.zeros(L)
+
+ # First condition
+ h_inv[0] = 1 / h[0]
+
+ # Recursive computation
+ for m in range(1, L):
+ s = 0
+ for i in range(1, m+1):
+ s += h_inv[m-i] * h[i]
+ h_inv[m] = -s / h[0]
+
+ return h_inv
+
+FIR_filter_coefficient = [[16, 20, 24, 28, 32, 40, 48],[1e-5, 1e-4, 1e-3, 1e-2, 1e-1],[1e-5, 1e-4, 1e-3, 1e-2, 1e-1]]
+inverse_FIR_filter_coefficient =[0.75,5e-2,'optimization' ]
+
+
+def analyze_and_plot_inverse_fir(
+ response: ArrayLike,
+ time: ArrayLike,
+ Ts: float = 0.5,
+ L_values: List[int] = FIR_filter_coefficient[0],
+ lam1_values: List[float] = FIR_filter_coefficient[1],
+ lam2_values: List[float] = FIR_filter_coefficient[2],
+ M: Optional[int] = None,
+ sigma_ns: float = inverse_FIR_filter_coefficient[0],
+ lam_smooth: float = inverse_FIR_filter_coefficient[1],
+ method: Literal['optimization', 'analytical'] = inverse_FIR_filter_coefficient[2],
+ verbose: bool = True
+ ) -> np.ndarray:
+ """
+ Analyze and plot the inverse FIR filter for a given FIR filter.
+
+ Parameters
+ ----------
+ response: ArrayLike
+ The response signal to analyze.
+ time: ArrayLike
+ The time array for the response signal.
+ Ts: float, optional
+ The sampling interval (in ns). Default is 0.5.
+ L_values: List[int], optional
+ The list of FIR filter lengths to search over. Default is [16, 20, 24, 28, 32, 40, 48].
+ lam1_values: List[float], optional
+ The list of regularization parameters for the identity matrix to search over. Default is [1e-5, 1e-4, 1e-3, 1e-2, 1e-1].
+ lam2_values: List[float], optional
+ The list of regularization parameters for the exponential tail to search over. Default is [1e-5, 1e-4, 1e-3, 1e-2, 1e-1].
+ M: int, optional
+ The length of the inverse FIR filter coefficients to solve for (default is len(h)).
+ sigma_ns: float, optional
+ The standard deviation of the target Gaussian for delta approximation (default is 0.75).
+ lam_smooth: float, optional
+ The regularization parameter for first-difference smoothing of h_inv (default is 5e-2).
+ method: Literal['optimization', 'analytical'], optional
+ The method to use for inversion ('optimization' or 'analytical'). Default is 'optimization'.
+ verbose: bool, optional
+ Whether to print verbose output. Default is True.
+
+ Returns
+ -------
+ best_h : ndarray
+ The best FIR filter coefficients.
+ h_inv : ndarray
+ The inverse FIR filter coefficients.
+ fig_fir_fit : matplotlib.figure.Figure
+ The figure object containing the FIR filter fit analysis and optimization plots.
+ fig : matplotlib.figure.Figure
+ The figure object containing the signal correction analysis and inverse FIR coefficients plots.
+ """
+ best_h, best_params, best_reconstructed, best_error, fig_fir_fit = analyze_and_plot_fir_fit(
+ response=response,
+ time=time,
+ Ts=Ts,
+ L_values=L_values,
+ lam1_values=lam1_values,
+ lam2_values=lam2_values,
+ verbose=verbose
+ )
+ h_inv = invert_fir(
+ h=best_h,
+ Ts=Ts,
+ M=M,
+ method=method,
+ sigma_ns=sigma_ns,
+ lam_smooth=lam_smooth
+ )
+ delta = conv_causal(best_h, h_inv, N=len(best_h))
+ ideal_response = np.ones(len(response))
+
+ # ===========================================
+ # Compute Corrected Signal
+ # ===========================================
+
+ # Method 1: Predistort ideal signal, then apply forward distortion
+ # (This simulates what would happen in practice)
+ L_guard = len(h_inv)
+ guard = np.zeros(L_guard)
+
+ # Pad ideal signal
+ ideal_padded = np.concatenate([guard, ideal_response, guard])
+
+ # Apply predistortion
+ predistorted_padded = conv_causal(ideal_padded, h_inv)
+
+ # Extract central region
+ start = L_guard
+ end = start + len(ideal_response)
+ predistorted_response = predistorted_padded[start:end]
+
+ # Apply forward distortion (simulate what happens in hardware)
+ corrected_response = conv_causal(predistorted_response, best_h, N=len(ideal_response))
+
+ # Compute correction error
+ correction_error = np.linalg.norm(corrected_response - ideal_response) / np.linalg.norm(ideal_response)
+ print(f"Correction error (NRMS): {correction_error:.3e}")
+
+ # Method 2: Apply inverse directly to measured signal
+ # (Alternative approach - corrects the measured signal directly)
+ distorted_padded = np.concatenate([guard, response, guard])
+ corrected_from_measured_padded = conv_causal(distorted_padded, h_inv)
+ corrected_from_measured = corrected_from_measured_padded[start:end]
+
+ correction_error_measured = np.linalg.norm(corrected_from_measured - ideal_response) / np.linalg.norm(ideal_response)
+ print(f"Correction error (from measured, NRMS): {correction_error_measured:.3e}")
+
+ # ===========================================
+ # Visualization
+ # ===========================================
+ print("\n" + "="*70)
+ print("Step 4: Visualization")
+ print("="*70)
+
+ fig, axes = plt.subplots(3, 2, figsize=(16, 12))
+
+ # Plot 1: Original signals
+ ax = axes[0, 0]
+ ax.plot(time, ideal_response, 'g--', label='Ideal Signal', linewidth=2, alpha=0.7)
+ ax.plot(time, response, 'r-', label='Distorted Signal', linewidth=2)
+ ax.plot(time, best_reconstructed, 'b:', label='Predicted (from FIR)', linewidth=2, alpha=0.7)
+ ax.axhline(y=1.001, color='gray', linestyle='--', linewidth=1, alpha=0.7)
+ ax.axhline(y=0.999, color='gray', linestyle='--', linewidth=1, alpha=0.7)
+ ax.set_ylim([0.95,1.05])
+ ax.set_xlabel('Time (ns)')
+ ax.set_ylabel('Amplitude')
+ ax.set_title('Original Signals and FIR Prediction')
+ ax.legend()
+ ax.grid(True, alpha=0.3)
+
+ # Plot 2: FIR Filters
+ ax = axes[0, 1]
+ ax.plot(best_h, 'b-o', label='Forward FIR (h)', markersize=4, linewidth=2)
+ ax.plot(h_inv, 'r-s', label='Inverse FIR (h_inv)', markersize=4, linewidth=2)
+ ax.set_xlabel('Tap Index')
+ ax.set_ylabel('Coefficient Value')
+ ax.set_title('Extracted FIR Filters')
+ ax.legend()
+ ax.grid(True, alpha=0.3)
+
+ # Plot 3: Predistorted and Corrected Signals
+ ax = axes[1, 0]
+ ax.plot(time, ideal_response, 'g--', label='Ideal Signal', linewidth=2, alpha=0.7)
+ ax.plot(time, predistorted_response, 'c-', label='Predistorted Signal', linewidth=2, alpha=0.7)
+ ax.plot(time, corrected_response, 'm-', label='Corrected Signal (sim)', linewidth=2)
+ ax.axhline(y=1.001, color='gray', linestyle='--', linewidth=1, alpha=0.7)
+ ax.axhline(y=0.999, color='gray', linestyle='--', linewidth=1, alpha=0.7)
+ ax.set_ylim([0.95,1.05])
+ ax.set_xlabel('Time (ns)')
+ ax.set_ylabel('Amplitude')
+ ax.set_title('Predistortion and Correction')
+ ax.legend()
+ ax.grid(True, alpha=0.3)
+
+ # Plot 4: Correction Comparison
+ ax = axes[1, 1]
+ ax.plot(time, ideal_response, 'g--', label='Ideal Signal', linewidth=2, alpha=0.7)
+ ax.plot(time, response, 'r-', label='Distorted Signal', linewidth=2, alpha=0.5)
+ ax.plot(time, corrected_response, 'm-', label='Corrected (predistort method)', linewidth=2)
+ ax.plot(time, corrected_from_measured, 'orange', linestyle='-',
+ label='Corrected (from measured)', linewidth=2, alpha=0.7)
+ ax.axhline(y=1.001, color='gray', linestyle='--', linewidth=1, alpha=0.7)
+ ax.axhline(y=0.999, color='gray', linestyle='--', linewidth=1, alpha=0.7)
+ ax.set_ylim([0.95,1.05])
+ ax.set_xlabel('Time (ns)')
+ ax.set_ylabel('Amplitude')
+ ax.set_title('Correction Comparison')
+ ax.legend()
+ ax.grid(True, alpha=0.3)
+
+ # Plot 5: Error Analysis
+ ax = axes[2, 0]
+ residual_fit = response - best_reconstructed
+ residual_correction = corrected_response - ideal_response
+ ax.plot(time, residual_fit, 'b-', label=f'Fit Residual (σ={np.std(residual_fit):.4e})', linewidth=1.5)
+ ax.plot(time, residual_correction, 'm-', label=f'Correction Residual (σ={np.std(residual_correction):.4e})', linewidth=1.5)
+ ax.axhline(y=0, color='k', linestyle='--', alpha=0.3)
+ ax.set_xlabel('Time (ns)')
+ ax.set_ylabel('Residual')
+ ax.set_title('Residual Analysis')
+ ax.legend()
+ ax.grid(True, alpha=0.3)
+
+ # Plot 6: h * h_inv (should be close to delta)
+ ax = axes[2, 1]
+ t_delta = np.arange(len(delta)) * Ts * 1e9 # Convert Ts (seconds) to ns for plotting
+ ax.plot(t_delta, delta, 'g-o', markersize=4, linewidth=2)
+ ax.axhline(y=0, color='k', linestyle='--', alpha=0.3)
+ ax.set_xlabel('Time (ns)')
+ ax.set_ylabel('Amplitude')
+ ax.set_title(f'h * h_inv (should approximate δ, peak={np.max(np.abs(delta)):.3e})')
+ ax.grid(True, alpha=0.3)
+
+ plt.tight_layout()
+ plt.show()
+ if verbose:
+ # ===========================================
+ # Summary
+ # ===========================================
+ print("\n" + "="*70)
+ print("SUMMARY")
+ print("="*70)
+ print(f"Forward FIR (h):")
+ print(f" - Length: {len(best_h)}")
+ print(f" - lam1: {best_params['lam1']:.0e}")
+ print(f" - lam2: {best_params['lam2']:.0e}")
+ print(f" - Sum: {np.sum(best_h):.6f}")
+ print(f" - Max coefficient: {np.max(np.abs(best_h)):.6f}")
+ print(f" - Fit error: {best_error:.3e} ({best_error*100:.2f}%)")
+
+ print(f"\nInverse FIR (h_inv):")
+ print(f" - Length: {len(h_inv)}")
+ print(f" - sigma_ns: {sigma_ns}")
+ print(f" - lam_smooth: {lam_smooth}")
+ print(f" - method: {method}")
+ print(f" - Sum: {np.sum(h_inv):.6f}")
+ print(f" - Max coefficient: {np.max(np.abs(h_inv)):.6f}")
+
+ print(f"\nCorrection Performance:")
+ print(f" - Correction error (predistort method): {correction_error:.3e} ({correction_error*100:.2f}%)")
+ print(f" - Correction error (from measured): {correction_error_measured:.3e} ({correction_error_measured*100:.2f}%)")
+
+ print(f"\nFilter Coefficients:")
+ print(f" Forward FIR (h): {best_h}..." if len(best_h) > 10 else f" Forward FIR (h): {best_h}")
+ print(f" Inverse FIR (h_inv): {h_inv}..." if len(h_inv) > 10 else f" Inverse FIR (h_inv): {h_inv}")
+ print("="*70)
+
+ return best_h, h_inv, fig_fir_fit, fig
+
+# %% example usage
+# load data
+with nc.Dataset('/Users/fabioansaloni/Downloads/ds (1).h5', 'r') as f:
+ t_data = np.array(f.variables['time'][:])
+ y_data = np.squeeze(np.array(f.variables['flux'][:]))
+
+# normalize the response
+normalized_response_raw = y_data / y_data[-10:].mean()
+# resample the response to 2 GS/s
+original_Ts = t_data[1] - t_data[0]
+normalized_response_2gsps = resample_to_target_rate(normalized_response_raw, original_Ts, 0.5)
+time_2gsps = np.arange(len(normalized_response_2gsps)) * 0.5 + 1
+# analyze and plot the forward and inverse FIR filters
+h_fir, inv_fir, fig_fir_fit, fig_inv_fir_fit = analyze_and_plot_inverse_fir(
+ response=normalized_response_2gsps,
+ time=time_2gsps
+ )
+
+# how to calculate the predistorted and corrected response
+ideal_response = np.ones(len(y_data))
+predistorted_response = lfilter(inv_fir, 1, ideal_response)
+corrected_response = lfilter(h_fir, 1, predistorted_response)
+
+# feedforward filters to be used in config file
+fir_filter = inv_fir.tolist() # or conv_causal(inv_fir, inv_fir_old).tolist() if inv_fir_old already exists
+# %%
diff --git a/qualang_tools/digital_filters/first_corr.png b/qualang_tools/digital_filters/first_corr.png
new file mode 100644
index 00000000..14ba772c
Binary files /dev/null and b/qualang_tools/digital_filters/first_corr.png differ
diff --git a/qualang_tools/digital_filters/secon_corr.png b/qualang_tools/digital_filters/secon_corr.png
new file mode 100644
index 00000000..312fb13c
Binary files /dev/null and b/qualang_tools/digital_filters/secon_corr.png differ
diff --git a/qualang_tools/digital_filters/step_IIR.png b/qualang_tools/digital_filters/step_IIR.png
new file mode 100644
index 00000000..7060c4f7
Binary files /dev/null and b/qualang_tools/digital_filters/step_IIR.png differ