How to calculate and use filter functions for arbitrary controls
Calculate the frequencydomain noise sensitivity of driven controls
Boulder Opal allows the calculation of filter functions as a way of estimating the sensitivity of an arbitrary control on a Ddimensional Hilbert space to timedependent noise expressed in the frequency domain. In this guide we show how to define, calculate, and visualize filter functions obtained using Boulder Opal.
For more general information on the mathematical underpinning of the filter function framework see Section III.B. of Software tools for quantum control: Improving quantum computer performance through noise and error suppression
Summary workflow
Usually, workflows involving filter function calculations include these four steps. A concrete example of these is in the next section.
1. Calculate filter functions
The function qctrl.functions.calculate_filter_function
takes the following parameters:

duration
, the duration of the control pulses, 
frequencies
, which is a list of frequencies where the filter function is to be sampled, 
drives
,shifts
, anddrifts
, which represent the different terms in the Hamiltonian (at least one term needs to be provided), 
sample_count
, which can be used to adjust the precision of the calculation of the filter functions (setting it toNone
uses the exact filter function calculation), 
projection_operator
, which projects into the subspace of interest, and 
result_scope
, which specifies whether to return the frequency domain noise operators.
More specifically, the three different types of Hamiltonian terms correspond to:

Drive
: a constant nonHermitian operator multiplied by a complex scalarvalued function of time, plus its Hermitian conjugate. It can be constructed using the objectqctrl.types.filter_function.Drive
. 
Shift
: a constant Hermitian operator multiplied by a real scalarvalued function of time that is constructed usingqctrl.types.filter_function.Shift
. 
Drift
: a constant Hermitian operator built usingqctrl.types.filter_function.Drift
.
2. Access the frequency domain noise operators
The frequency domain noise operators are computed as part of the filter function calculation. They are available in the result samples when the result_scope
parameter is set appropriately.
3. Visualize filter functions
After their calculation, the numerical values of the filter functions are stored at filter_function_result.samples
(where filter_function_result
is the object returned by qctrl.functions.calculate_filter_function
). Each sample contains a corresponding frequency
(our x coordinates), and an inverse_power
(our y coordinates).
This information can be plotted using the plot_filter_functions
function from the QCTRL Visualizer Python package.
The resulting graph represents the “noise admittance” of the control—where it is low, the noise is suppressed and where it is high, the noise is transmitted.
4. Calculate control infidelity due to noise with filter functions
As explained in the reference documentation for qctrl.functions.calculate_filter_function
, the overlap integrals of the filter functions with the noise power spectral densities (PSDs) can be used to estimate the operational infidelity in the weak noise regime.
If multiple sources of noise are present, the total infidelity $\mathcal{I}$ can be estimated by adding the contributions of each noise as represented by each filter function $F_k(f)$ and power spectral density $S_k(f)$ where $k$ is the index for the list of filter functions:
With the filter functions we calculated previously, we can find the operational infidelity with respect to power spectral densities that we define. To do this, we just need to numerically calculate the overlap integrals, which can be done using the np.trapz
NumPy function for composite trapezoidal rule integration. In general this firstorder approximation provides a good quantitative estimate for noiseinfidelities $\lesssim 10\%$. At higher noise strengths simple extensions may be used to approximate higherorder contributions (see Experimental noise filtering by quantum control, Supplementary information).
Note that if the pulses were not capable of perfectly creating the target X gates in the absence of noise, we should add the value of the noisefree infidelity, $\mathcal{I}_0$, to the calculation of the total infidelity: $ \mathcal{I} \approx \mathcal{I}_0 + \sum_k \int_{\infty}^\infty F_k (f) S_k (f) \, \mathrm{d}f$. In the case of pulses that were predefined to generate the desired gate, $\mathcal{I}_0 \approx 0$.
Worked example: Noise sensitivity of composite pulses applied to a single qubit
In this example, we will use filter functions to compare the sensitivity of different composite $\pi$pulses to timevarying amplitude and dephasing noise. The Hamiltonian of the system we will be considering is:
\[H(t) = \frac{1 + \beta_\Omega(t)}{2} \left[ \Omega(t) \sigma_ + \Omega^* (t) \sigma_+ \right] + \frac{\Delta(t)}{2} \sigma_z + \frac{\eta(t)}{2} \sigma_z,\]where $\Omega(t)$ is a timedependent Rabi rate, $\beta_\Omega(t)$ is a fractional timedependent amplitude fluctuation process, $\Delta(t)$ is a timedependent clock shift, $\eta(t)$ is a small, slowlyvarying stochastic dephasing noise process, and $\sigma_k$ are the Pauli matrices (with $\sigma_\pm = (\sigma_x \mp i \sigma_y)/2$).
Because there are two different noise channels we will be calculating two different filter functions for each control we study. This will permit a direct comparison of how each control exhibits different responses to different noise channels. We will compare the filter functions of each for a range of noise frequencies from $10^{8} \Omega_\mathrm{max}$ to $\Omega_\mathrm{max}$, where $\Omega_\mathrm{max}/2\pi = 1~\mathrm{MHz}$ is the maximum Rabi frequency.
Create the pulses
We will consider the following driven control schemes for the controllable $\Omega (t)$ and $\Delta (t)$ terms:
 Primitive: sensitive to both dephasing and amplitude noise
 BB1: suppresses amplitude noise
 CORPSE: suppresses dephasing noise
 CORPSE in BB1: suppresses both dephasing and amplitude noise
These schemes are available from QCTRL Open Controls and described in the reference documentation.
We first set up Python objects representing the pulse, control, and noise. In the Hamiltonian used in this example, the drive corresponds to the term that contains the Rabi rate ($\Omega(t)$ and $\Omega^*(t)$), the shift corresponds to the term that contains the clock shift ($\Delta(t)$), and the drift corresponds to the dephasing term.
The Rabi coupling term and the clockshift term are defined using a pulse from QCTRL Open Controls.
As we are interested in separately studying the robustness against amplitude noise in the complexvalued controls and in the dephasing term, we define two versions of the drive: one with noise, and one without it.
import matplotlib.pyplot as plt
import numpy as np
from attr import asdict
from dataclasses import dataclass
from typing import List
# Predefined pulse imports
from qctrlopencontrols import (
new_bb1_control,
new_corpse_control,
new_corpse_in_bb1_control,
new_primitive_control,
)
from qctrlopencontrols.driven_controls.driven_control import DrivenControl
from qctrlvisualizer import plot_filter_functions
from qctrl import Qctrl
# Starting a session with the API
qctrl = Qctrl()
# Define standard matrices
identity = np.array([[1, 0], [0, 1]], dtype=complex)
sigma_z = np.array([[1, 0], [0, 1]], dtype=complex)
sigma_x = np.array([[0, 1], [1, 0]], dtype=complex)
sigma_m = np.array([[0, 1], [0, 0]], dtype=complex)
# Define control parameters
omega_max = 2 * np.pi * 1e6 # Hz
total_rotation = np.pi
# Define schemes for driven controls to compare
pulses: List[DrivenControl] = [
pulse_function(
rabi_rotation=total_rotation,
azimuthal_angle=0.0,
maximum_rabi_rate=omega_max,
name=name,
)
for name, pulse_function in [
("primitive", new_primitive_control),
("BB1", new_bb1_control),
("CORPSE", new_corpse_control),
("CORPSE in BB1", new_corpse_in_bb1_control),
]
]
@dataclass
class ControlScheme:
name: str
amplitude_filter_function_samples: qctrl.types.filter_function.Sample
dephasing_filter_function_samples: qctrl.types.filter_function.Sample
def create_control_scheme(
pulse: DrivenControl, frequencies: np.ndarray
) > ControlScheme:
# Define noiseless controls
noiseless_drive = qctrl.types.filter_function.Drive(
control=[
qctrl.types.ComplexSegmentInput(duration=duration, value=value)
for duration, value in zip(
pulse.durations, pulse.rabi_rates * np.exp(1j * pulse.azimuthal_angles)
)
],
operator=sigma_m / 2,
)
noiseless_shift = qctrl.types.filter_function.Shift(
control=[
qctrl.types.RealSegmentInput(duration=duration, value=value)
for duration, value in zip(pulse.durations, pulse.detunings)
],
operator=sigma_z / 2,
)
# Define complex control with multiplicative noise
noisy_drive = qctrl.types.filter_function.Drive(
control=noiseless_drive.control, operator=noiseless_drive.operator, noise=True
)
# Define additive noise
dephasing_noise = qctrl.types.filter_function.Drift(
noise=True, operator=sigma_z / 2
)
# The filter function of the amplitude noise excludes the term
# for the dephasing noise (the drift)
amplitude_result = qctrl.functions.calculate_filter_function(
duration=pulse.duration,
frequencies=frequencies,
sample_count=3000,
drives=[noisy_drive],
shifts=[noiseless_shift],
result_scope="NO_FREQUENCY_DOMAIN_NOISE_OPERATORS",
)
assert not amplitude_result.errors
# The filter function of the dephasing noise excludes the term
# for the amplitude noise (the noisy drive)
dephasing_result = qctrl.functions.calculate_filter_function(
duration=pulse.duration,
frequencies=frequencies,
sample_count=None,
drives=[noiseless_drive],
shifts=[noiseless_shift],
drifts=[dephasing_noise],
result_scope="ALL",
)
assert not dephasing_result.errors
return ControlScheme(
name=pulse.name,
amplitude_filter_function_samples=amplitude_result.samples,
dephasing_filter_function_samples=dephasing_result.samples,
)
Calculate and visualize the filter functions
Here we calculate filter functions and set the frequencies in the array to be spaced logarithmically, because we will also be interested in plotting the filter functions on a loglog graph.
# Define filter function parameters
interpolated_frequencies = omega_max * np.logspace(8, 0, 1000, base=10)
# We expect two calls to calculate_filter_function for each pulse as we are evaluating
# for both drifts and shifts
control_schemes = [
create_control_scheme(pulse, interpolated_frequencies) for pulse in pulses
]
Your task calculate_filter_function (action_id="1164537") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_filter_function (action_id="1164537") has completed.
Your task calculate_filter_function (action_id="1164553") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_filter_function (action_id="1164553") has completed.
Your task calculate_filter_function (action_id="1164558") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_filter_function (action_id="1164558") has completed.
Your task calculate_filter_function (action_id="1164564") has completed.
Your task calculate_filter_function (action_id="1164570") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_filter_function (action_id="1164570") has completed.
Your task calculate_filter_function (action_id="1164576") has completed.
Your task calculate_filter_function (action_id="1164580") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_filter_function (action_id="1164580") has completed.
Your task calculate_filter_function (action_id="1164584") has completed.
Amplitude noise
plot_filter_functions(
plt.figure(),
{
scheme.name: map(asdict, scheme.amplitude_filter_function_samples)
for scheme in control_schemes
},
)
We see that the BB1 and CORPSEinBB1 pulses show small values of the filter function at low frequencies, corresponding to their known amplitudenoisesuppressing characteristics. CORPSE and primitive gates do not exhibit any suppression of amplitude noise.
Dephasing noise
plot_filter_functions(
plt.figure(),
{
scheme.name: map(asdict, scheme.dephasing_filter_function_samples)
for scheme in control_schemes
},
)
We see that the CORPSE and CORPSEinBB1 pulses show small values of the filter function at low frequencies, corresponding to their known dephasingnoisesuppressing characteristics. BB1 and primitive gates do not exhibit any suppression of dephasing noise.
Calculate the infidelity with respect to noise power spectral densities using noise operators
The frequency domain noise operators are computed as part of the filter function calculation. They are available in the result samples when the result_scope
parameter is set appropriately. In this case we have configured the filter function calculation for the dephasing noise to include the operators.
With these in hand we can calculate the expected infidelity of a control in the presence of noise channels exhibiting specific quantitative power spectra. We define dimensionful power spectra below and compare the performance of the pulses.
scheme = [s for s in control_schemes if s.name == "primitive"][0]
last_dephasing_sample = scheme.dephasing_filter_function_samples[1]
print(f"Scheme: {scheme.name}")
print(f"Frequency domain noise operator at f={last_dephasing_sample.frequency}Hz")
print(last_dephasing_sample.frequency_domain_noise_operator)
Scheme: primitive
Frequency domain noise operator at f=6283185.307179586Hz
[[1.00946617e082.11765689e08j 1.60661531e09+3.37035562e09j]
[1.60661531e093.37035562e09j 1.00946617e08+2.11765689e08j]]
# Define amplitude noise power spectral density
psd_amplitude = lambda frequency: 1e2 / (frequency**2 + 1)
# Define dephasing noise power spectral density
psd_dephasing = lambda frequency: 1e12 / (frequency**2 + 1)
for scheme in control_schemes:
# Calculate overlap integral for amplitude noise
frequencies = np.array(
[sample.frequency for sample in scheme.amplitude_filter_function_samples]
)
inverse_powers = np.array(
[sample.inverse_power for sample in scheme.amplitude_filter_function_samples]
)
amplitude_infidelity = np.trapz(
y=inverse_powers * psd_amplitude(frequencies), x=frequencies
)
# Calculate overlap integral for dephasing noise
frequencies = np.array(
[sample.frequency for sample in scheme.dephasing_filter_function_samples]
)
inverse_powers = np.array(
[sample.inverse_power for sample in scheme.dephasing_filter_function_samples]
)
dephasing_infidelity = np.trapz(
y=inverse_powers * psd_dephasing(frequencies), x=frequencies
)
# Add contributions from the two noises
infidelity = amplitude_infidelity + dephasing_infidelity
# Print results
print()
print(f"Scheme: {scheme.name}")
print(f"Infidelity due to amplitude noise: {amplitude_infidelity:.2e}")
print(f"Infidelity due to dephasing noise: {dephasing_infidelity:.2e}")
print(f"Total infidelity: {infidelity:.2e}")
Scheme: primitive
Infidelity due to amplitude noise: 3.72e02
Infidelity due to dephasing noise: 3.82e02
Total infidelity: 7.54e02
Scheme: BB1
Infidelity due to amplitude noise: 7.50e07
Infidelity due to dephasing noise: 3.82e02
Total infidelity: 3.82e02
Scheme: CORPSE
Infidelity due to amplitude noise: 3.73e02
Infidelity due to dephasing noise: 7.50e07
Total infidelity: 3.73e02
Scheme: CORPSE in BB1
Infidelity due to amplitude noise: 2.84e06
Infidelity due to dephasing noise: 2.19e06
Total infidelity: 5.02e06
Summary
The plots and infidelities show that the BB1 controls perform better against (lowfrequency) amplitude noise than the primitive $\pi$pulse, but perform just like the primitive in the case of dephasing noise. This is expected, as BB1 is one of the controlerrorcompensating driven controls. The inverse is true of the pure CORPSE controls: they perform as poorly as the primitive against amplitude noise, but perform better against (lowfrequency) dephasing noise. This behavior is expected as CORPSE is a dephasingerrorcompensating driven control. Finally, the dephasingandcontrolerrorcompensating driven control CORPSE in BB1 performs better than the primitive for lowfrequencies of both noises (albeit less so than the controls specialized for one specific kind of noise).
We have thus demonstrated how Boulder Opal can be used to characterize the sensitivity of different controls to timedependent noise channels by calculating their corresponding filter functions.
This notebook was run using the following package versions. It should also be compatible with newer versions of the QCTRL Python package.
Package  version 

Python  3.9.12 
matplotlib  3.5.1 
numpy  1.21.5 
scipy  1.7.3 
qctrl  19.1.0 
qctrlcommons  17.1.1 
qctrlopencontrols  9.1.0 
qctrltoolkit  1.5.0 
qctrlvisualizer  3.2.1 