# Filter functions

**Calculation of noise filtering properties in driven controls**

The Q-CTRL Python package allows the calculation of filter functions as a way of estimating the sensitivity of a control to the frequency of a time-dependent noise. In this guide we show how to define, calculate, and visualize filter functions obtained using the Q-CTRL Python package.

```
# Essential imports
from attr import asdict
import numpy as np
from qctrl import Qctrl
# Predefined pulse imports
from qctrlopencontrols import new_predefined_driven_control
# Plotting imports
import matplotlib.pyplot as plt
from qctrlvisualizer import plot_filter_functions
# Starting a session with the API
qctrl = Qctrl()
```

## Worked example: composite $\pi$ -pulses applied to a single qubit under amplitude and dephasing noise

In this example, we will compare the filter functions corresponding to different composite $\pi$-pulses under 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 time-dependent Rabi rate, $\beta_\Omega(t)$ is a fractional time-dependent amplitude fluctuation process, $\Delta(t)$ is a time-dependent clock shift, $\eta(t)$ is a small, slowly-varying stochastic dephasing noise process, and $\sigma_k$ are the Pauli matrices (with $\sigma_\pm = (\sigma_x \mp i \sigma_y)/2$).

We will consider the following driven control schemes for the controllable $\Omega (t)$ and $\Delta (t)$ terms:

- primitive,
- BB1,
- CORPSE,
- CORPSE in BB1.

These schemes are available from Q-CTRL Open Controls and described in the reference documentation.

In this system, filter functions can be calculated for two kinds of noise, so that the sensitivity of the controls against them can be compared. 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.

### Creating the pulses

As described in the Setup user guide, we first set up Python objects representing the pulse, control, and noise.

The `qctrl.functions.calculate_filter_function`

function accepts the Hamiltonian in the form of separate terms for each of the complex-valued controls (`drives`

), real-valued controls (`shifts`

), and static terms (`drifts`

). 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 clock shift term are defined using a pulse from Q-CTRL Open Controls.

As we are interested in studying separately the robustness against amplitude noise in the complex-valued controls and in the dephasing term, we define two versions of the drive: one with noise, and one without it.

```
# Define standard matrices
identity = np.array([[1, 0], [0, 1]], dtype=np.complex)
sigma_z = np.array([[1, 0], [0, -1]], dtype=np.complex)
sigma_x = np.array([[0, 1], [1, 0]], dtype=np.complex)
sigma_m = np.array([[0, 1], [0, 0]], dtype=np.complex)
# Define control parameters
omega_max = 2 * np.pi * 1e6 # Hz
total_rotation = np.pi
# Define schemes for driven controls to compare
schemes = {scheme: {} for scheme in ['primitive', 'BB1', 'CORPSE', 'CORPSE in BB1']}
for scheme, scheme_objects in schemes.items():
# Define pulse objects using pulses from Q-CTRL Open Controls
pulse = new_predefined_driven_control(
rabi_rotation=total_rotation,
azimuthal_angle=0.,
maximum_rabi_rate=omega_max,
scheme=scheme,
name=scheme,
)
# 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,
)
# Save relevant objects for later use
scheme_objects['total_duration'] = pulse.duration
scheme_objects['noiseless_drives'] = [noiseless_drive]
scheme_objects['noiseless_shifts'] = [noiseless_shift]
scheme_objects['noisy_drives'] = [noisy_drive]
scheme_objects['noisy_drifts'] = [dephasing_noise]
```

### Calculating the filter functions

Filter functions are created using `qctrl.functions.calculate_filter_function`

, which 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`

, and`drifts`

, 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,`projection_operator`

, which projects into the subspace of interest, and`result_scope`

, which specifies whether to return the frequency domain noise operators .

In this example, the frequencies in the array will be spaced logarithmically, because we will also be interested in plotting the filter functions on a log-log graph.

```
# Define filter function parameters
sample_count = 3000
interpolated_frequencies = omega_max*np.logspace(-8, 0, 1000, base=10)
# Create filter function objects
for scheme_objects in schemes.values():
# The filter function of the amplitude noise excludes the term
# for the dephasing noise (the drift).
scheme_objects['amplitude_filter_function'] = qctrl.functions.calculate_filter_function(
duration=scheme_objects['total_duration'],
frequencies=interpolated_frequencies,
sample_count=sample_count,
drives=scheme_objects['noisy_drives'],
shifts=scheme_objects['noiseless_shifts'],
result_scope="NO_FREQUENCY_DOMAIN_NOISE_OPERATORS",
)
# The filter function of the dephasing noise excludes the term
# for the amplitude noise (the noisy drive).
scheme_objects['dephasing_filter_function'] = qctrl.functions.calculate_filter_function(
duration=scheme_objects['total_duration'],
frequencies=interpolated_frequencies,
sample_count=sample_count,
drives=scheme_objects['noiseless_drives'],
shifts=scheme_objects['noiseless_shifts'],
drifts=scheme_objects['noisy_drifts'],
result_scope="ALL",
)
```

### Visualizing the 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 Q-CTRL Python Visualizer package.

```
plot_filter_functions(
plt.figure(),
{
scheme: [asdict(sample) for sample in scheme_objects['amplitude_filter_function'].samples]
for scheme, scheme_objects in schemes.items()
}
)
```

```
plot_filter_functions(
plt.figure(),
{
scheme: [asdict(sample) for sample in scheme_objects['dephasing_filter_function'].samples]
for scheme, scheme_objects in schemes.items()
}
)
```

### Accessing 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. In this case we have configured the filter function calculation for the dephasing noise to include the operators.

```
print("Scheme: {}".format('primitive'))
print("Frequency domain noise operator at f =",schemes['primitive']['dephasing_filter_function'].samples[-1].frequency," Hz:")
print(schemes['primitive']['dephasing_filter_function'].samples[-1].frequency_domain_noise_operator)
```

### Calculating the infidelity with respect to noise power spectral densities

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)\}$:

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 calculate numerically the overlap integrals, which can be done using the `np.trapz`

NumPy function for composite trapezoidal rule integration.

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 noise-free 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$. For other cases, you can calculate their noise-free infidelity by following the instructions from the Simulation user guide.

```
# Define amplitude noise power spectral density
psd_amplitude = lambda frequency: 1e-2/(frequency**2 + 1)
# Define dephasing noise power spectral density
psd_dephasing = lambda frequency: 1e12/(frequency**2 + 1)
for scheme, scheme_objects in schemes.items():
# Calculate overlap integral for amplitude noise
frequencies = np.array([
sample.frequency
for sample in scheme_objects['amplitude_filter_function'].samples
])
inverse_powers = np.array([
sample.inverse_power
for sample in scheme_objects['amplitude_filter_function'].samples
])
scheme_objects["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_objects['dephasing_filter_function'].samples
])
inverse_powers = np.array([
sample.inverse_power
for sample in scheme_objects['dephasing_filter_function'].samples
])
scheme_objects["dephasing_infidelity"] = np.trapz(
y=inverse_powers*psd_dephasing(frequencies),
x=frequencies,
)
# Add contributions from the two noises
scheme_objects["infidelity"] = (
scheme_objects["amplitude_infidelity"] + scheme_objects["dephasing_infidelity"]
)
# Print results
print("\nScheme: {}".format(scheme))
print("Infidelity due to amplitude noise: {}".format(scheme_objects["amplitude_infidelity"]))
print("Infidelity due to dephasing noise: {}".format(scheme_objects["dephasing_infidelity"]))
print("Total infidelity: {}".format(scheme_objects["infidelity"]))
```

### Summary

The plots and infidelities show that the BB1 controls perform better against (low-frequency) 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 *control-error-compensating 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 (low-frequency) dephasing noise. This behavior is expected as CORPSE is a *dephasing-error-compensating driven control*. Finally, the *dephasing-and-control-error-compensating driven control* CORPSE in BB1 performs better than the primitive for low-frequencies of both noises (albeit less so than the controls specialized for one specific kind of noise).

We have thus demonstrated how the Q-CTRL Python package can be used to characterize the sensitivity of different controls to time-dependent noise channels by calculating their corresponding filter functions.