diff --git a/filter_functions/basis.py b/filter_functions/basis.py
index 08f2da9..506be2b 100644
--- a/filter_functions/basis.py
+++ b/filter_functions/basis.py
@@ -568,7 +568,7 @@ def _full_from_partial(elems: Sequence, traceless: bool, labels: Sequence[str])
     else:
         ggm = Basis.ggm(elems.d)
 
-    coeffs = expand(elems, ggm, tidyup=True)
+    coeffs = expand(elems, ggm, hermitian=elems.isherm, tidyup=True)
     # Throw out coefficient vectors that are all zero (should only happen for
     # the identity)
     coeffs = coeffs[(coeffs != 0).any(axis=-1)]
@@ -637,7 +637,7 @@ def normalize(b: Basis) -> Basis:
 
 
 def expand(M: Union[ndarray, Basis], basis: Union[ndarray, Basis],
-           normalized: bool = True, tidyup: bool = False) -> ndarray:
+           normalized: bool = True, hermitian: bool = False, tidyup: bool = False) -> ndarray:
     r"""
     Expand the array *M* in the basis given by *basis*.
 
@@ -650,6 +650,9 @@ def expand(M: Union[ndarray, Basis], basis: Union[ndarray, Basis],
         The basis of shape (m, d, d) in which to expand.
     normalized: bool {True}
         Wether the basis is normalized.
+    hermitian: bool (default: False)
+        If M is hermitian along its last two axes, the result will be
+        real.
     tidyup: bool {False}
         Whether to set values below the floating point eps to zero.
 
@@ -673,15 +676,18 @@ def expand(M: Union[ndarray, Basis], basis: Union[ndarray, Basis],
                     {\mathrm{tr}\big(C_j^\dagger C_j\big)}.
 
     """
-    coefficients = np.tensordot(M, basis, axes=[(-2, -1), (-1, -2)])
 
+    def cast(arr): return arr.real if hermitian and basis.isherm else arr
+
+    coefficients = cast(np.tensordot(M, basis, axes=[(-2, -1), (-1, -2)]))
     if not normalized:
-        coefficients /= np.einsum('bij,bji->b', basis, basis).real
+        coefficients /= cast(np.einsum('bij,bji->b', basis, basis))
 
     return util.remove_float_errors(coefficients) if tidyup else coefficients
 
 
-def ggm_expand(M: Union[ndarray, Basis], traceless: bool = False) -> ndarray:
+def ggm_expand(M: Union[ndarray, Basis], traceless: bool = False,
+               hermitian: bool = False) -> ndarray:
     r"""
     Expand the matrix *M* in a Generalized Gell-Mann basis [Bert08]_.
     This function makes use of the explicit construction prescription of
@@ -699,11 +705,14 @@ def ggm_expand(M: Union[ndarray, Basis], traceless: bool = False) -> ndarray:
         expansion. If it is known beforehand that M is traceless, the
         corresponding coefficient is zero and thus doesn't need to be
         computed.
+    hermitian: bool (default: False)
+        If M is hermitian along its last two axes, the result will be
+        real.
 
     Returns
     -------
     coefficients: ndarray
-        The coefficient array with shape (d**2,) or (..., d**2)
+        The real coefficient array with shape (d**2,) or (..., d**2)
 
     References
     ----------
@@ -715,6 +724,8 @@ def ggm_expand(M: Union[ndarray, Basis], traceless: bool = False) -> ndarray:
     if M.shape[-1] != M.shape[-2]:
         raise ValueError('M should be square in its last two axes')
 
+    def cast(arr): return arr.real if hermitian else arr
+
     # Squeeze out an extra dimension to be shape agnostic
     square = M.ndim < 3
     if square:
@@ -740,18 +751,19 @@ def ggm_expand(M: Union[ndarray, Basis], traceless: bool = False) -> ndarray:
     diag_idx_shifted = tuple([...] + [tuple(i for i in range(1, d))]*2)
 
     # Compute the coefficients
-    coeffs = np.zeros((*M.shape[:-2], d**2), dtype=complex)
+    coeffs = np.zeros((*M.shape[:-2], d**2), dtype=float if hermitian else complex)
     if not traceless:
         # First element is proportional to the trace of M
-        coeffs[..., 0] = np.einsum('...jj', M)/np.sqrt(d)
+        coeffs[..., 0] = cast(np.einsum('...jj', M))/np.sqrt(d)
 
     # Elements proportional to the symmetric GGMs
-    coeffs[..., sym_rng] = (M[triu_idx] + M[tril_idx])/np.sqrt(2)
+    coeffs[..., sym_rng] = cast(M[triu_idx] + M[tril_idx])/np.sqrt(2)
     # Elements proportional to the antisymmetric GGMs
-    coeffs[..., sym_rng + n_sym] = 1j*(M[triu_idx] - M[tril_idx])/np.sqrt(2)
+    coeffs[..., sym_rng + n_sym] = cast(1j*(M[triu_idx] - M[tril_idx]))/np.sqrt(2)
     # Elements proportional to the diagonal GGMs
-    coeffs[..., diag_rng + 2*n_sym] = M[diag_idx].cumsum(axis=-1) - diag_rng*M[diag_idx_shifted]
-    coeffs[..., diag_rng + 2*n_sym] /= np.sqrt(diag_rng*(diag_rng + 1))
+    coeffs[..., diag_rng + 2*n_sym] = cast(M[diag_idx].cumsum(axis=-1)
+                                           - diag_rng*M[diag_idx_shifted])
+    coeffs[..., diag_rng + 2*n_sym] /= cast(np.sqrt(diag_rng*(diag_rng + 1)))
 
     return coeffs.squeeze() if square else coeffs
 
diff --git a/filter_functions/gradient.py b/filter_functions/gradient.py
index 0817f4a..1ab9eb8 100644
--- a/filter_functions/gradient.py
+++ b/filter_functions/gradient.py
@@ -174,11 +174,12 @@ def _liouville_derivative(dt: Coefficients, propagators: ndarray, basis: Basis,
 
     # Calculate the whole propagator derivative
     propagators_deriv = np.zeros((c_opers_transformed.shape[1], n-1, n, d, d), dtype=complex)
+    U_deriv_transformed = np.zeros((c_opers_transformed.shape[1], n-1, d, d), dtype=complex)
     for g in range(n-1):
-        propagators_deriv[:, g, :g+1] = (propagators[g+1]
-                                         @ propagators[1:g+2].conj().swapaxes(-1, -2)
-                                         @ U_deriv[:, :g+1]
-                                         @ propagators[:g+1])
+        U_deriv_transformed[:, g] = (propagators[g+1].conj().swapaxes(-1, -2)
+                                     @ U_deriv[:, g]
+                                     @ propagators[g])
+        propagators_deriv[:, g, :g+1] = propagators[g+1] @ U_deriv_transformed[:, :g+1]
 
     # Equivalent but usually slower contraction: 'htsba,jbc,tcd,kda->thsjk'
     # Can just take 2*Re(·) when calculating x + x*
@@ -466,8 +467,8 @@ def calculate_derivative_of_control_matrix_from_scratch(
 
     See Also
     --------
-    :func:`_liouville_derivative`
-    :func:`_control_matrix_at_timestep_derivative`
+    _liouville_derivative
+    _control_matrix_at_timestep_derivative
     """
     # Distinction between control and drift operators and only
     # calculate the derivatives in control direction
diff --git a/filter_functions/numeric.py b/filter_functions/numeric.py
index aad714b..68c8408 100644
--- a/filter_functions/numeric.py
+++ b/filter_functions/numeric.py
@@ -444,8 +444,8 @@ def calculate_noise_operators_from_atomic(
 
     See Also
     --------
-    :func:`calculate_noise_operators_from_scratch`: Compute the operators from scratch.
-    :func:`calculate_control_matrix_from_atomic`: Same calculation in Liouville space.
+    calculate_noise_operators_from_scratch: Compute the operators from scratch.
+    calculate_control_matrix_from_atomic: Same calculation in Liouville space.
     """
     n = len(noise_operators_atomic)
     # Allocate memory
@@ -572,8 +572,8 @@ def calculate_noise_operators_from_scratch(
 
     See Also
     --------
-    :func:`calculate_noise_operators_from_atomic`: Compute the operators from atomic segments.
-    :func:`calculate_control_matrix_from_scratch`: Same calculation in Liouville space.
+    calculate_noise_operators_from_atomic: Compute the operators from atomic segments.
+    calculate_control_matrix_from_scratch: Same calculation in Liouville space.
     """
     if t is None:
         t = np.concatenate(([0], np.asarray(dt).cumsum()))
@@ -624,7 +624,7 @@ def calculate_noise_operators_from_scratch(
     return noise_operators
 
 
-@util.parse_optional_parameters({'which': ('total', 'correlations')})
+@util.parse_optional_parameters(which=('total', 'correlations'))
 def calculate_control_matrix_from_atomic(
         phases: ndarray,
         control_matrix_atomic: ndarray,
@@ -859,7 +859,7 @@ def calculate_control_matrix_from_scratch(
 
 def calculate_control_matrix_periodic(phases: ndarray, control_matrix: ndarray,
                                       total_propagator_liouville: ndarray,
-                                      repeats: int) -> ndarray:
+                                      repeats: int, check_invertible: bool = True) -> ndarray:
     r"""
     Calculate the control matrix of a periodic pulse given the phase
     factors, control matrix and transfer matrix of the total propagator,
@@ -877,6 +877,11 @@ def calculate_control_matrix_periodic(phases: ndarray, control_matrix: ndarray,
         pulse.
     repeats: int
         The number of repetitions.
+    check_invertible: bool, optional
+        Check for all frequencies if the inverse :math:`\mathbb{I} -
+        e^{i\omega T} \mathcal{Q}^{(1)}` exists by calculating the
+        determinant. If it does not exist, the sum is evaluated
+        explicitly for those points. The default is True.
 
     Returns
     -------
@@ -909,28 +914,23 @@ def calculate_control_matrix_periodic(phases: ndarray, control_matrix: ndarray,
     # evaluate it explicitly.
     eye = np.eye(total_propagator_liouville.shape[0])
     T = np.multiply.outer(phases, total_propagator_liouville)
-
-    # Mask the invertible frequencies. The chosen atol is somewhat empiric.
     M = eye - T
-    M_inv = nla.inv(M)
-    good_inverse = np.isclose(M_inv @ M, eye, atol=1e-10, rtol=0).all((1, 2))
+    if check_invertible:
+        invertible = ~np.isclose(nla.det(M), 0)
+    else:
+        invertible = np.array(True)
 
     S = np.empty((*phases.shape, *total_propagator_liouville.shape), dtype=complex)
-    # Evaluate the sum for invertible frequencies
-    S[good_inverse] = M_inv[good_inverse] @ (eye - nla.matrix_power(T[good_inverse], repeats))
-
-    # Evaluate the sum for non-invertible frequencies
-    # HINT: Using numba, this could be a factor of ten or so faster. But since
-    # usually very few omega-values have a bad inverse, the compilation
-    # overhead is not compensated.
-    if (~good_inverse).any():
-        S[~good_inverse] = eye + sum(accumulate(repeat(T[~good_inverse], repeats-1), np.matmul))
+    # Solve LSE instead of computing inverse, faster + numerically more stable
+    S[invertible] = nla.solve(M[invertible], eye - nla.matrix_power(T[invertible], repeats))
+    if (~invertible).any():
+        S[~invertible] = eye + sum(accumulate(repeat(T[~invertible], repeats-1), np.matmul))
 
     control_matrix_tot = (control_matrix.transpose(2, 0, 1) @ S).transpose(1, 2, 0)
     return control_matrix_tot
 
 
-@util.parse_optional_parameters({'which': ('total', 'correlations')})
+@util.parse_optional_parameters(which=('total', 'correlations'))
 def calculate_cumulant_function(
         pulse: 'PulseSequence',
         spectrum: Optional[ndarray] = None,
@@ -1096,23 +1096,25 @@ def calculate_cumulant_function(
     if d == 2 and pulse.basis.btype in ('Pauli', 'GGM'):
         # Single qubit case. Can use simplified expression
         cumulant_function = np.zeros(decay_amplitudes.shape, decay_amplitudes.dtype)
-        diag_mask = np.eye(N, dtype=bool)
+        diag_mask = np.zeros((N, N), dtype=bool)
+        diag_mask[1:, 1:] = ~np.eye(N-1, dtype=bool)
 
         # Offdiagonal terms
-        cumulant_function[..., ~diag_mask] = decay_amplitudes[..., ~diag_mask]
+        cumulant_function[..., diag_mask] = decay_amplitudes[..., diag_mask]
 
         # Diagonal terms K_ii given by sum over diagonal of Gamma excluding
         # Gamma_ii. Since the Pauli basis is traceless, K_00 is zero, therefore
         # start at K_11.
-        diag_idx = deque((True, False, True, True))
+        diag_deque = deque((False, True, True))
         for i in range(1, N):
+            diag_idx = [False] + list(diag_deque)
             cumulant_function[..., i, i] = -decay_amplitudes[..., diag_idx, diag_idx].sum(axis=-1)
             # shift the item not summed over by one
-            diag_idx.rotate()
+            diag_deque.rotate()
 
         if second_order:
-            cumulant_function -= frequency_shifts
-            cumulant_function += frequency_shifts.swapaxes(-1, -2)
+            cumulant_function[..., 1:, 1:] -= frequency_shifts[..., 1:, 1:]
+            cumulant_function[..., 1:, 1:] += frequency_shifts[..., 1:, 1:].swapaxes(-1, -2)
     else:
         # Multi qubit case. Use general expression. Drop imaginary part since
         # result is guaranteed to be real (if we didn't do anything wrong)
@@ -1134,7 +1136,7 @@ def calculate_cumulant_function(
     return cumulant_function.real
 
 
-@util.parse_optional_parameters({'which': ('total', 'correlations')})
+@util.parse_optional_parameters(which=('total', 'correlations'))
 def calculate_decay_amplitudes(
         pulse: 'PulseSequence',
         spectrum: ndarray,
@@ -1353,7 +1355,7 @@ def calculate_frequency_shifts(
     return frequency_shifts
 
 
-@util.parse_which_FF_parameter
+@util.parse_optional_parameters(which=('fidelity', 'generalized'))
 def calculate_filter_function(control_matrix: ndarray, which: str = 'fidelity') -> ndarray:
     r"""Compute the filter function from the control matrix.
 
@@ -1524,8 +1526,10 @@ def calculate_second_order_filter_function(
         intermediates = dict()
 
     # Work around possibly populated intermediates dict with missing keys
-    n_opers_transformed = intermediates.get('n_opers_transformed',
-                                            _transform_hamiltonian(eigvecs, n_opers, n_coeffs))
+    n_opers_transformed = intermediates.get('n_opers_transformed')
+    if n_opers_transformed is None:
+        n_opers_transformed = _transform_hamiltonian(eigvecs, n_opers, n_coeffs)
+
     try:
         basis_transformed_cache = intermediates['basis_transformed']
         ctrlmat_step_cache = intermediates['control_matrix_step']
@@ -1585,7 +1589,7 @@ def calculate_second_order_filter_function(
     return result
 
 
-@util.parse_which_FF_parameter
+@util.parse_optional_parameters(which=('fidelity', 'generalized'))
 def calculate_pulse_correlation_filter_function(control_matrix: ndarray,
                                                 which: str = 'fidelity') -> ndarray:
     r"""Compute pulse correlation filter function from control matrix.
@@ -1786,7 +1790,7 @@ def error_transfer_matrix(
     For non-Gaussian noise the expression above is perturbative and
     includes noise up to order :math:`\xi^2` and hence
     :math:`\tilde{\mathcal{U}} = \mathbb{1} + \mathcal{K}(\tau) +
-    \mathcal{O}(\xi^2)`
+    \mathcal{O}(\xi^4)`
     (although it is evaluated as a matrix exponential in any case).
 
     Given the above expression of the error transfer matrix, the
@@ -1829,7 +1833,7 @@ def error_transfer_matrix(
     return error_transfer_matrix
 
 
-@util.parse_optional_parameters({'which': ('total', 'correlations')})
+@util.parse_optional_parameters(which=('total', 'correlations'))
 def infidelity(
         pulse: 'PulseSequence',
         spectrum: Union[Coefficients, Callable],
@@ -2028,12 +2032,12 @@ def infidelity(
         else:
             raise ValueError("spacing should be either 'linear' or 'log'.")
 
-        delta_n = (n_max - n_min)//n_points
+        delta_n = (n_max - n_min)//(n_points - 1)
         n_samples = np.arange(n_min, n_max + delta_n, delta_n)
 
         convergence_infids = np.empty((len(n_samples), len(idx)))
         for i, n in enumerate(n_samples):
-            freqs = xspace(omega_IR, omega_UV, n//2)
+            freqs = xspace(omega_IR, omega_UV, n)
             convergence_infids[i] = infidelity(pulse, spectrum(freqs), freqs,
                                                n_oper_identifiers=n_oper_identifiers,
                                                which='total', show_progressbar=show_progressbar,
diff --git a/filter_functions/plotting.py b/filter_functions/plotting.py
index 786bb98..f390d1c 100644
--- a/filter_functions/plotting.py
+++ b/filter_functions/plotting.py
@@ -85,14 +85,17 @@ def _make_str_tex_compatible(s: str) -> str:
     return s
 
 
-def get_bloch_vector(states: Sequence[State]) -> ndarray:
-    r"""
-    Get the Bloch vector from quantum states.
-    """
+def _import_qutip_or_raise():
     try:
         import qutip as qt
     except ImportError as err:
         raise RuntimeError('Requirements not fulfilled. Please install Qutip') from err
+    return qt
+
+
+def get_bloch_vector(states: Sequence[State]) -> ndarray:
+    """Get the Bloch vector from quantum states."""
+    qt = _import_qutip_or_raise()
 
     if isinstance(states[0], qt.Qobj):
         a = np.empty((3, len(states)))
@@ -110,10 +113,7 @@ def get_bloch_vector(states: Sequence[State]) -> ndarray:
 def init_bloch_sphere(**bloch_kwargs) -> qt.Bloch:
     """A helper function to create a Bloch instance with a default viewing
     angle and axis labels."""
-    try:
-        import qutip as qt
-    except ImportError as err:
-        raise RuntimeError('Requirements not fulfilled. Please install Qutip') from err
+    qt = _import_qutip_or_raise()
 
     bloch_kwargs.setdefault('view', [-150, 30])
     b = qt.Bloch(**bloch_kwargs)
@@ -127,7 +127,7 @@ def init_bloch_sphere(**bloch_kwargs) -> qt.Bloch:
     return b
 
 
-@util.parse_optional_parameters({'prop': ['total', 'piecewise']})
+@util.parse_optional_parameters(prop=('total', 'piecewise'))
 def get_states_from_prop(U: Sequence[Operator], psi0: Optional[State] = None,
                          prop: str = 'total') -> ndarray:
     r"""
@@ -228,10 +228,9 @@ def plot_bloch_vector_evolution(
         b = init_bloch_sphere(fig=fig, axes=axes, **bloch_kwargs)
 
     if n_samples is None:
-        # 5 time points during  the smallest time interval in pulse.t. Being
-        # careful that doesn't blow up in our face for extremely narrow pulses,
-        # max out at 5000.
-        n_samples = min([5000, 5*int(pulse.tau/np.diff(pulse.t).min())])
+        # At least 100, at most 5000 points, default 10 points per smallest
+        # time interval
+        n_samples = min(5000, max(10*int(pulse.tau/pulse.dt.min()), 100))
 
     times = np.linspace(pulse.t[0], pulse.tau, n_samples)
     propagators = pulse.propagator_at_arb_t(times)
@@ -624,7 +623,8 @@ def plot_pulse_correlation_filter_function(
     return fig, axes, legend
 
 
-def plot_infidelity_convergence(n_samples: Sequence[int], infids: Sequence[float]) -> FigureAxes:
+def plot_infidelity_convergence(n_samples: Sequence[int], infids: Sequence[float],
+                                axes: Optional[Axes] = None) -> FigureAxes:
     """
     Plot the convergence of the infidelity integral. The function
     arguments are those returned by
@@ -639,6 +639,8 @@ def plot_infidelity_convergence(n_samples: Sequence[int], infids: Sequence[float
     infids: array_like, shape (n_samples, [n_oper_inds, optional])
         Array with the calculated infidelities for each noise operator
         on the second axis or the second axis already traced out.
+    axes: sequence of two matplotlib axes, optional
+        Two axes that the result is plotted in.
 
     Returns
     -------
@@ -648,21 +650,25 @@ def plot_infidelity_convergence(n_samples: Sequence[int], infids: Sequence[float
         The matplotlib axes instances used for plotting.
 
     """
-    fig, ax = plt.subplots(2, 1, sharex=True)
-    ax[1].set_xlabel(r'$n_\omega$')
-    ax[0].set_ylabel(r'$\mathcal{I}$')
-    ax[1].set_ylabel(r'$|\Delta\mathcal{I}|/\mathcal{I}$ (%)')
-    ax[0].set_xlim(min(n_samples), max(n_samples))
-    ax[0].grid()
-    ax[1].grid()
+    if axes is None:
+        fig, axes = plt.subplots(2, 1, sharex=True)
+    else:
+        fig = axes[0].get_figure()
+
+    axes[1].set_xlabel(r'$n_\omega$')
+    axes[0].set_ylabel(r'$\mathcal{I}$')
+    axes[1].set_ylabel(r'$|\Delta\mathcal{I}|/\mathcal{I}$ (%)')
+    axes[0].set_xlim(min(n_samples), max(n_samples))
+    axes[0].grid()
+    axes[1].grid()
 
-    ax[0].plot(n_samples, infids, 'o-')
-    ax[1].semilogy(n_samples, np.abs(np.gradient(infids, axis=0))/infids*100, 'o-')
+    axes[0].plot(n_samples, infids, 'o-')
+    axes[1].semilogy(n_samples, np.abs(np.gradient(infids, axis=0))/infids*100, 'o-')
 
-    return fig, ax
+    return fig, axes
 
 
-@util.parse_optional_parameters({'colorscale': ['linear', 'log']})
+@util.parse_optional_parameters(colorscale=('linear', 'log'))
 def plot_cumulant_function(
         pulse: Optional['PulseSequence'] = None,
         spectrum: Optional[ndarray] = None,
@@ -795,7 +801,7 @@ def plot_cumulant_function(
         else:
             basis_labels = [f'$C_{{{i}}}$' for i in range(K.shape[-1])]
     else:
-        if basis_labels and len(basis_labels) != K.shape[-1]:
+        if basis_labels is not None and len(basis_labels) != K.shape[-1]:
             raise ValueError('Invalid number of basis_labels given')
 
         basis_labels = [_make_str_tex_compatible(bl) for bl in basis_labels]
diff --git a/filter_functions/pulse_sequence.py b/filter_functions/pulse_sequence.py
index 2082c24..b376311 100644
--- a/filter_functions/pulse_sequence.py
+++ b/filter_functions/pulse_sequence.py
@@ -40,7 +40,7 @@
 """
 
 import bisect
-from copy import copy
+import copy
 from itertools import accumulate, compress, zip_longest
 from typing import Any, Dict, Iterable, List, Mapping, Optional, Sequence, Tuple, Union
 from warnings import warn
@@ -287,7 +287,7 @@ def __init__(self, *args, **kwargs) -> None:
 
     def __str__(self):
         """String method."""
-        return f'{repr(self)} with total duration {self.tau}'
+        return f'{repr(self)}\n\tof dimension {self.d} and duration {self.duration}'
 
     def __eq__(self, other: object) -> bool:
         """
@@ -393,9 +393,15 @@ def __getitem__(self, key) -> 'PulseSequence':
     def __copy__(self) -> 'PulseSequence':
         """Return shallow copy of self"""
         cls = self.__class__
-        copy = cls.__new__(cls)
-        copy.__dict__.update(self.__dict__)
-        return copy
+        copied = cls.__new__(cls)
+        copied.__dict__.update(self.__dict__)
+        return copied
+
+    def __deepcopy__(self, memo=None) -> 'PulseSequence':
+        cls = self.__class__
+        copied = cls.__new__(cls)
+        copied.__dict__.update({key: copy.deepcopy(val) for key, val in self.__dict__.items()})
+        return copied
 
     def __matmul__(self, other: 'PulseSequence') -> 'PulseSequence':
         """Concatenation of PulseSequences."""
@@ -586,7 +592,7 @@ def get_pulse_correlation_control_matrix(self) -> ndarray:
             "True."
         )
 
-    @util.parse_optional_parameters({'which': ('fidelity', 'generalized'), 'order': (1, 2)})
+    @util.parse_optional_parameters(which=('fidelity', 'generalized'), order=(1, 2))
     def get_filter_function(
             self, omega: Coefficients,
             which: str = 'fidelity',
@@ -686,7 +692,7 @@ def get_filter_function(
             # order == 2
             return self._filter_function_2
 
-    @util.parse_optional_parameters({'which': ('fidelity', 'generalized'), 'order': (1, 2)})
+    @util.parse_optional_parameters(which=('fidelity', 'generalized'), order=(1, 2))
     def cache_filter_function(
             self, omega: Coefficients,
             control_matrix: Optional[ndarray] = None,
@@ -748,6 +754,7 @@ def cache_filter_function(
                         self._filter_function_pc = F_pc
                     else:
                         # which == 'generalized'
+                        self._filter_function_pc = F_pc.trace(axis1=4, axis2=5)
                         self._filter_function_pc_gen = F_pc
 
                     filter_function = F_pc.sum(axis=(0, 1))
@@ -767,12 +774,13 @@ def cache_filter_function(
                 self._filter_function = filter_function
             else:
                 # which == 'generalized'
+                self._filter_function = filter_function.trace(axis1=2, axis2=3)
                 self._filter_function_gen = filter_function
         else:
             # order == 2
             self._filter_function_2 = filter_function
 
-    @util.parse_which_FF_parameter
+    @util.parse_optional_parameters(which=('fidelity', 'generalized'))
     def get_pulse_correlation_filter_function(self, which: str = 'fidelity') -> ndarray:
         r"""
         Get the pulse correlation filter function given by
@@ -1005,8 +1013,8 @@ def nbytes(self) -> int:
 
         return sum(_nbytes)
 
-    @util.parse_optional_parameters({'method': ('conservative', 'greedy',
-                                                'frequency dependent', 'all')})
+    @util.parse_optional_parameters(method=('conservative', 'greedy',
+                                            'frequency dependent', 'all'))
     def cleanup(self, method: str = 'conservative') -> None:
         """
         Delete cached byproducts of the calculation of the filter
@@ -1306,8 +1314,8 @@ def _concatenate_Hamiltonian(
 
     if any(len(value) > 1 for value in oper_to_identifier_mapping.values()):
         # Clash: two different identifiers are assigned to the same operator
-        raise ValueError('Trying to concatenate pulses with equal operators with different ' +
-                         'identifiers. Please choose unique identifiers!')
+        raise ValueError(f'Trying to concatenate pulses with equal {kind} operators but ' +
+                         f'different identifiers. Please choose unique {kind} identifiers!')
 
     # A dict that maps the identifiers of each Hamiltonian to the identifiers
     # in the new Hamiltonian
@@ -1561,7 +1569,7 @@ def concatenate_without_filter_function(pulses: Iterable[PulseSequence],
     return newpulse
 
 
-@util.parse_which_FF_parameter
+@util.parse_optional_parameters(which=('fidelity', 'generalized'))
 def concatenate(
         pulses: Iterable[PulseSequence],
         calc_pulse_correlation_FF: bool = False,
@@ -1620,7 +1628,7 @@ def concatenate(
     """
     pulses = tuple(pulses)
     if len(pulses) == 1:
-        return copy(pulses[0])
+        return copy.deepcopy(pulses[0])
 
     newpulse, _, n_oper_mapping = concatenate_without_filter_function(
         pulses, return_identifier_mappings=True
@@ -1804,8 +1812,6 @@ def concatenate_periodic(pulse: PulseSequence, repeats: int) -> PulseSequence:
     except TypeError:
         raise TypeError(f'Expected pulses to be iterable, not {type(pulse)}')
 
-    cached_ctrl_mat = pulse.is_cached('control_matrix')
-
     # Initialize a new PulseSequence instance with the Hamiltonians sequenced
     # (this is much easier than in the general case, thus do it on the fly)
     dt = np.tile(pulse.dt, repeats)
@@ -1822,7 +1828,7 @@ def concatenate_periodic(pulse: PulseSequence, repeats: int) -> PulseSequence:
     )
     newpulse.tau = repeats*pulse.tau
 
-    if not cached_ctrl_mat:
+    if not pulse.is_cached('control_matrix'):
         # No cached filter functions to reuse and pulse correlation FFs not
         # requested. If they were, continue even if there are no cached FF
         # they cannot be computed anymore afterwards.
@@ -1830,16 +1836,13 @@ def concatenate_periodic(pulse: PulseSequence, repeats: int) -> PulseSequence:
 
     phases_at = pulse.get_total_phases(pulse.omega)
     control_matrix_at = pulse.get_control_matrix(pulse.omega)
-    L_at = pulse.total_propagator_liouville
+    total_propagator_liouville_at = pulse.total_propagator_liouville
 
     newpulse.total_propagator = nla.matrix_power(pulse.total_propagator, repeats)
     newpulse.cache_total_phases(pulse.omega)
-    # Might be cheaper for small repeats to use matrix_power, but this function
-    # is aimed at a large number so we calculate it explicitly
-    newpulse.total_propagator_liouville = newpulse.total_propagator_liouville
 
     control_matrix_tot = numeric.calculate_control_matrix_periodic(phases_at, control_matrix_at,
-                                                                   L_at, repeats)
+                                                                   total_propagator_liouville_at, repeats)
 
     newpulse.cache_filter_function(pulse.omega, control_matrix_tot)
 
diff --git a/filter_functions/superoperator.py b/filter_functions/superoperator.py
index 9422f05..4e755d7 100644
--- a/filter_functions/superoperator.py
+++ b/filter_functions/superoperator.py
@@ -77,19 +77,13 @@ def liouville_representation(U: ndarray, basis: _b.Basis) -> ndarray:
     Hilbert space.
     """
     U = np.asanyarray(U)
+    conjugated_basis = np.einsum('...ba,ibc,...cd->...iad', U.conj(), basis, U,
+                                 optimize=['einsum_path', (1, 2), (0, 1)])
     if basis.btype == 'GGM' and basis.d > 12:
         # Can do closed form expansion and overhead compensated
-        path = ['einsum_path', (0, 1), (0, 1)]
-        conjugated_basis = np.einsum('...ba,ibc,...cd->...iad', U.conj(), basis, U, optimize=path)
-        # If the basis is hermitian, the result will be strictly real so we can
-        # drop the imaginary part
-        R = _b.ggm_expand(conjugated_basis).real
+        return _b.ggm_expand(conjugated_basis, hermitian=basis.isherm)
     else:
-        path = ['einsum_path', (0, 1), (0, 1), (0, 1)]
-        R = np.einsum('...ba,ibc,...cd,jda', U.conj(), basis, U, basis,
-                      optimize=path).real
-
-    return R
+        return _b.expand(conjugated_basis, basis, hermitian=basis.isherm)
 
 
 def liouville_to_choi(superoperator: ndarray, basis: _b.Basis) -> ndarray:
diff --git a/filter_functions/util.py b/filter_functions/util.py
index 3dab060..2a7175a 100644
--- a/filter_functions/util.py
+++ b/filter_functions/util.py
@@ -75,7 +75,6 @@
 from typing import Callable, Dict, Iterable, List, Optional, Sequence, Tuple, Union
 
 import numpy as np
-from numpy import linalg as nla
 from numpy import ndarray
 
 from .types import Operator, State
@@ -153,18 +152,20 @@ def cexp(x: ndarray, out=None, where=True) -> ndarray:
     return out
 
 
-def parse_optional_parameters(params_dict: Dict[str, Sequence]) -> Callable:
-    """Decorator factory to parse optional parameter with certain legal values.
+def parse_optional_parameters(**allowed_kwargs: Dict[str, Sequence]) -> Callable:
+    """Decorator factory to parse optional parameter with certain legal
+    values.
 
-    For ``params_dict = {name: allowed, ...}``: If the parameter value
-    corresponding to ``name`` (either in args or kwargs of the decorated
-    function) is not contained in ``allowed`` a ``ValueError`` is raised.
+    For ``allowed_kwargs = {name: allowed, ...}``: If the parameter
+    value corresponding to ``name`` (either in args or kwargs of the
+    decorated function) is not contained in ``allowed`` a ``ValueError``
+    is raised.
     """
     def decorator(func):
         @functools.wraps(func)
         def wrapper(*args, **kwargs):
             parameters = inspect.signature(func).parameters
-            for name, allowed in params_dict.items():
+            for name, allowed in allowed_kwargs.items():
                 idx = tuple(parameters).index(name)
                 try:
                     value = args[idx]
@@ -180,9 +181,6 @@ def wrapper(*args, **kwargs):
     return decorator
 
 
-parse_which_FF_parameter = parse_optional_parameters({'which': ('fidelity', 'generalized')})
-
-
 def parse_operators(opers: Sequence[Operator], err_loc: str) -> List[ndarray]:
     """Parse a sequence of operators and convert to ndarray.
 
@@ -902,7 +900,6 @@ def oper_equiv(psi: Union[Operator, State], phi: Union[Operator, State],
     """
     # Convert qutip.Qobj's to numpy arrays
     psi, phi = [obj.full() if hasattr(obj, 'full') else obj for obj in (psi, phi)]
-
     psi, phi = np.atleast_2d(psi, phi)
 
     if eps is None:
@@ -915,14 +912,15 @@ def oper_equiv(psi: Union[Operator, State], phi: Union[Operator, State],
             eps *= (np.prod(psi.shape[-2:])*phi.shape[-1]*2)**2
 
     try:
-        inner_product = (psi.swapaxes(-1, -2).conj() @ phi).trace(0, -1, -2)
+        # Don't need to round at this point
+        inner_product = dot_HS(psi, phi, eps=0)
     except ValueError as err:
         raise ValueError('psi and phi have incompatible dimensions!') from err
 
     if normalized:
         norm = 1
     else:
-        norm = nla.norm(psi, axis=(-1, -2))*nla.norm(phi, axis=(-1, -2))
+        norm = np.sqrt(dot_HS(psi, psi, eps=0)*dot_HS(phi, phi, eps=0))
 
     phase = np.angle(inner_product)
     modulus = abs(inner_product)
@@ -940,6 +938,9 @@ def dot_HS(U: Operator, V: Operator, eps: Optional[float] = None) -> float:
     ----------
     U, V: qutip.Qobj or ndarray
         Objects to compute the inner product of.
+    eps: float
+        The floating point precision. The result is rounded to
+        `abs(int(np.log10(eps)))` decimals if `eps > 0`.
 
     Returns
     -------
@@ -970,17 +971,15 @@ def dot_HS(U: Operator, V: Operator, eps: Optional[float] = None) -> float:
             eps = 0
 
     if eps == 0:
-        decimals = 0
+        res = np.einsum('...ij,...ij', U.conj(), V)
     else:
-        decimals = abs(int(np.log10(eps)))
+        res = np.around(np.einsum('...ij,...ij', U.conj(), V), decimals=abs(int(np.log10(eps))))
 
-    res = np.round(np.einsum('...ij,...ij', U.conj(), V), decimals)
     return res if res.imag.any() else res.real
 
 
-@parse_optional_parameters({'spacing': ('log', 'linear')})
-def get_sample_frequencies(pulse: 'PulseSequence', n_samples: int = 300,
-                           spacing: str = 'log',
+@parse_optional_parameters(spacing=('log', 'linear'))
+def get_sample_frequencies(pulse: 'PulseSequence', n_samples: int = 300, spacing: str = 'log',
                            include_quasistatic: bool = False) -> ndarray:
     """Get *n_samples* sample frequencies spaced 'linear' or 'log'.
 
diff --git a/tests/test_basis.py b/tests/test_basis.py
index 575c3b3..5fc0e90 100644
--- a/tests/test_basis.py
+++ b/tests/test_basis.py
@@ -158,6 +158,21 @@ def test_basis_properties(self):
 
     def test_basis_expansion_and_normalization(self):
         """Correct expansion of operators and normalization of bases"""
+        # dtype
+        b = ff.Basis.ggm(3)
+        r = ff.basis.expand(rng.standard_normal((3, 3)), b, hermitian=False)
+        self.assertTrue(r.dtype == 'complex128')
+        r = ff.basis.expand(testutil.rand_herm(3), b, hermitian=True)
+        self.assertTrue(r.dtype == 'float64')
+        b._isherm = False
+        r = ff.basis.expand(testutil.rand_herm(3), b, hermitian=True)
+        self.assertTrue(r.dtype == 'complex128')
+
+        r = ff.basis.ggm_expand(testutil.rand_herm(3), hermitian=True)
+        self.assertTrue(r.dtype == 'float64')
+        r = ff.basis.ggm_expand(rng.standard_normal((3, 3)), hermitian=False)
+        self.assertTrue(r.dtype == 'complex128')
+
         for _ in range(10):
             d = rng.integers(2, 16)
             ggm_basis = ff.Basis.ggm(d)
diff --git a/tests/test_core.py b/tests/test_core.py
index f05add4..245ce85 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -21,8 +21,8 @@
 """
 This module tests the core functionality of the package.
 """
+import copy
 import string
-from copy import copy
 from random import sample
 
 import numpy as np
@@ -217,9 +217,6 @@ def test_pulse_sequence_constructor(self):
         pulse
         print(pulse)
 
-        # Hit __copy__ method
-        _ = copy(pulse)
-
         # Fewer identifiers than opers
         pulse_2 = ff.PulseSequence(
             [[util.paulis[1], [1], 'X'],
@@ -231,6 +228,24 @@ def test_pulse_sequence_constructor(self):
         self.assertArrayEqual(pulse_2.c_oper_identifiers, ('A_1', 'X'))
         self.assertArrayEqual(pulse_2.n_oper_identifiers, ('B_0', 'Y'))
 
+    def test_copy(self):
+        pulse = testutil.rand_pulse_sequence(2, 2)
+        old_copers = pulse.c_opers.copy()
+
+        copied = copy.copy(pulse)
+        deepcopied = copy.deepcopy(pulse)
+
+        self.assertEqual(pulse, copied)
+        self.assertEqual(pulse, deepcopied)
+
+        pulse.c_opers[...] = rng.standard_normal(size=pulse.c_opers.shape)
+
+        self.assertArrayEqual(pulse.c_opers, copied.c_opers)
+        self.assertArrayEqual(old_copers, deepcopied.c_opers)
+
+        self.assertEqual(pulse, copied)
+        self.assertNotEqual(pulse, deepcopied)
+
     def test_pulse_sequence_attributes(self):
         """Test attributes of single instance"""
         X, Y, Z = util.paulis[1:]
@@ -613,6 +628,21 @@ def test_cache_intermediates(self):
         self.assertArrayAlmostEqual(pulse._intermediates['basis_transformed'], basis_transformed,
                                     atol=1e-14)
 
+    def test_cache_filter_function(self):
+        omega = rng.random(32)
+        pulse = testutil.rand_pulse_sequence(2, 3, n_nops=2)
+        F_fidelity = numeric.calculate_filter_function(pulse.get_control_matrix(omega),
+                                                       'fidelity')
+        F_generalized = numeric.calculate_filter_function(pulse.get_control_matrix(omega),
+                                                          'generalized')
+
+        pulse.cache_filter_function(omega, filter_function=F_generalized, which='generalized')
+        self.assertTrue(pulse.is_cached('filter function'))
+        self.assertTrue(pulse.is_cached('generalized filter function'))
+
+        self.assertArrayEqual(pulse.get_filter_function(omega, which='generalized'), F_generalized)
+        self.assertArrayEqual(pulse.get_filter_function(omega, which='fidelity'), F_fidelity)
+
     def test_filter_function(self):
         """Test the filter function calculation and related methods"""
         for d, n_dt in zip(rng.integers(2, 10, (3,)),
@@ -776,6 +806,8 @@ def test_pulse_correlation_filter_function(self):
         pulse_3 = ff.concatenate([pulses['X'], pulses['Y']],
                                  calc_pulse_correlation_FF=True,
                                  which='generalized')
+        pulse_4 = copy.copy(pulse_3)
+        pulse_4.cleanup('all')
 
         self.assertTrue(pulse_2.is_cached('control_matrix_pc'))
         self.assertTrue(pulse_2.is_cached('filter_function_pc'))
@@ -829,6 +861,36 @@ def test_pulse_correlation_filter_function(self):
             )
         )
 
+        # Test caching
+        pulse_4.cache_filter_function(omega, control_matrix=control_matrix_pc, which='fidelity')
+        self.assertTrue(pulse_4.is_cached('pulse correlation control matrix'))
+        self.assertTrue(pulse_4.is_cached('pulse correlation filter function'))
+        self.assertTrue(pulse_4.is_cached('filter function'))
+        self.assertArrayAlmostEqual(pulse_3.get_pulse_correlation_control_matrix(),
+                                    pulse_4.get_pulse_correlation_control_matrix())
+        self.assertArrayAlmostEqual(pulse_3.get_pulse_correlation_filter_function(),
+                                    pulse_4.get_pulse_correlation_filter_function())
+        self.assertArrayAlmostEqual(pulse_3.get_filter_function(omega),
+                                    pulse_4.get_filter_function(omega))
+
+        pulse_4.cleanup('all')
+        pulse_4.cache_filter_function(omega, control_matrix=control_matrix_pc, which='generalized')
+        self.assertTrue(pulse_4.is_cached('pulse correlation control matrix'))
+        self.assertTrue(pulse_4.is_cached('generalized pulse correlation filter function'))
+        self.assertTrue(pulse_4.is_cached('generalized filter function'))
+        self.assertTrue(pulse_4.is_cached('pulse correlation filter function'))
+        self.assertTrue(pulse_4.is_cached('filter function'))
+        self.assertArrayAlmostEqual(pulse_3.get_pulse_correlation_control_matrix(),
+                                    pulse_4.get_pulse_correlation_control_matrix())
+        self.assertArrayAlmostEqual(pulse_3.get_pulse_correlation_filter_function('fidelity'),
+                                    pulse_4.get_pulse_correlation_filter_function('fidelity'))
+        self.assertArrayAlmostEqual(pulse_3.get_pulse_correlation_filter_function('generalized'),
+                                    pulse_4.get_pulse_correlation_filter_function('generalized'))
+        self.assertArrayAlmostEqual(pulse_3.get_filter_function(omega, which='fidelity'),
+                                    pulse_4.get_filter_function(omega, which='fidelity'))
+        self.assertArrayAlmostEqual(pulse_3.get_filter_function(omega, which='generalized'),
+                                    pulse_4.get_filter_function(omega, which='generalized'))
+
         # If for some reason filter_function_pc_xy is removed, check if
         # recovered from control_matrix_pc
         pulse_2._filter_function_pc = None
diff --git a/tests/test_plotting.py b/tests/test_plotting.py
index dd70591..3a9260a 100644
--- a/tests/test_plotting.py
+++ b/tests/test_plotting.py
@@ -260,7 +260,7 @@ def test_plot_cumulant_function(self):
         K = numeric.calculate_cumulant_function(complicated_pulse, spectrum, omega)
         fig, grid = plotting.plot_cumulant_function(
             complicated_pulse, spectrum=spectrum, omega=omega,
-            n_oper_identifiers=n_oper_identifiers, basis_labels=basis_labels,
+            n_oper_identifiers=n_oper_identifiers, basis_labels=np.array(basis_labels),
             basis_labelsize=4, linthresh=1e-4, cmap=plt.cm.jet
         )
         fig, grid = plotting.plot_cumulant_function(
@@ -306,13 +306,13 @@ def test_plot_cumulant_function(self):
         plt.close('all')
 
     def test_plot_infidelity_convergence(self):
-        def spectrum(omega):
-            return omega**0
-
-        n, infids = ff.infidelity(simple_pulse, spectrum, {},
-                                  test_convergence=True)
+        n, infids = ff.infidelity(simple_pulse, lambda x: x**0, {}, test_convergence=True)
         fig, ax = plotting.plot_infidelity_convergence(n, infids)
 
+        fig, ax = plt.subplots(1, 2)
+        cfig, ax = plotting.plot_infidelity_convergence(n, infids, ax)
+        self.assertIs(cfig, fig)
+
 
 class LaTeXRenderingTest(testutil.TestCase):
 
diff --git a/tests/test_precision.py b/tests/test_precision.py
index d55163e..0def1ea 100644
--- a/tests/test_precision.py
+++ b/tests/test_precision.py
@@ -549,6 +549,11 @@ def test_single_qubit_error_transfer_matrix(self):
         d = 2
         for n_dt in rng.integers(1, 11, 10):
             pulse = testutil.rand_pulse_sequence(d, n_dt, 3, 2, btype='Pauli')
+            traceless = rng.integers(2, dtype=bool)
+            if not traceless:
+                # Test that result correct for finite trace n_oper (edge case)
+                pulse.n_opers[0] = np.eye(d)/np.sqrt(d)
+
             omega = util.get_sample_frequencies(pulse, n_samples=51)
             n_oper_identifiers = pulse.n_oper_identifiers
             traces = pulse.basis.four_element_traces.todense()
@@ -575,14 +580,16 @@ def test_single_qubit_error_transfer_matrix(self):
                       np.einsum('...kl,kijl->...ij', Gamma, traces))/2
                 U_onfoot = sla.expm(K.sum(tuple(range(K.ndim - 2))))
                 U_from_K = ff.error_transfer_matrix(cumulant_function=K)
-                I_fidelity = ff.infidelity(pulse, S, omega)
-                I_decayamps = -np.einsum('...ii', K)/d**2
-                I_transfer = 1 - np.einsum('...ii', U)/d**2
                 self.assertArrayAlmostEqual(Up, U)
-                self.assertArrayAlmostEqual(I_fidelity, I_decayamps)
-                self.assertArrayAlmostEqual(I_transfer, I_fidelity.sum(), rtol=1e-4, atol=1e-10)
                 self.assertArrayAlmostEqual(U, U_onfoot, atol=1e-14)
                 self.assertArrayAlmostEqual(U_from_K, U_onfoot)
+                if traceless:
+                    # Simplified fidelity calculation relies on traceless n_opers
+                    I_fidelity = ff.infidelity(pulse, S, omega)
+                    I_decayamps = -np.einsum('...ii', K)/d**2
+                    I_transfer = 1 - np.einsum('...ii', U)/d**2
+                    self.assertArrayAlmostEqual(I_fidelity, I_decayamps)
+                    self.assertArrayAlmostEqual(I_transfer, I_fidelity.sum(), rtol=1e-4, atol=1e-10)
 
                 # second order
                 K -= (np.einsum('...kl,klji->...ij', Delta, traces) -
@@ -592,10 +599,13 @@ def test_single_qubit_error_transfer_matrix(self):
                 U = ff.error_transfer_matrix(pulse, S, omega, second_order=True)
                 U_onfoot = sla.expm(K.sum(tuple(range(K.ndim - 2))))
                 U_from_K = ff.error_transfer_matrix(cumulant_function=K)
-                I_transfer = 1 - np.einsum('...ii', U)/d**2
-                self.assertArrayAlmostEqual(I_transfer, I_fidelity.sum(), rtol=1e-4, atol=1e-10)
                 self.assertArrayAlmostEqual(U, U_onfoot, atol=1e-14)
                 self.assertArrayAlmostEqual(U_from_K, U_onfoot)
+                if traceless:
+                    # Simplified fidelity calculation relies on traceless n_opers
+                    I_transfer = 1 - np.einsum('...ii', U)/d**2
+                    self.assertArrayAlmostEqual(I_transfer, I_fidelity.sum(), rtol=1e-4, atol=1e-10)
+
 
     def test_multi_qubit_error_transfer_matrix(self):
         """Test the calculation of the multi-qubit transfer matrix"""
diff --git a/tests/test_sequencing.py b/tests/test_sequencing.py
index 0e30411..1032492 100644
--- a/tests/test_sequencing.py
+++ b/tests/test_sequencing.py
@@ -24,7 +24,7 @@
 
 import string
 from copy import deepcopy
-from itertools import product
+from itertools import product, repeat
 from random import sample
 
 import numpy as np
@@ -507,7 +507,7 @@ def test_different_n_opers(self):
             # Test forcibly caching
             self.assertTrue(pulse_13_2.is_cached('filter_function'))
 
-    def test_concatenation_periodic(self):
+    def test_concatenate_periodic(self):
         """Test concatenation for periodic Hamiltonians"""
         X, Y, Z = util.paulis[1:]
         A = 0.01
@@ -567,6 +567,24 @@ def test_concatenation_periodic(self):
         self.assertArrayAlmostEqual(F_LAB, F_CC, atol=1e-13)
         self.assertArrayAlmostEqual(F_LAB, F_CC_PERIODIC, atol=1e-13)
 
+        # Random tests, make sure we test for G=1
+        for d, G in zip(rng.integers(2, 7, 11), np.concatenate([[1], rng.integers(2, 1000, 10)])):
+            pulse = testutil.rand_pulse_sequence(d, 5, 2, 2)
+            pulse.cache_filter_function(rng.random(37))
+            a = ff.concatenate(repeat(pulse, G))
+            b = ff.concatenate_periodic(pulse, G)
+            self.assertEqual(a, b)
+            self.assertArrayAlmostEqual(a._control_matrix, b._control_matrix)
+            self.assertArrayAlmostEqual(a._filter_function, b._filter_function)
+
+            cm = ff.numeric.calculate_control_matrix_periodic(pulse.get_total_phases(pulse.omega),
+                                                              pulse.get_control_matrix(pulse.omega),
+                                                              pulse.total_propagator_liouville,
+                                                              G, check_invertible=False)
+            # Check mostly always equal
+            self.assertGreater(np.isclose(cm, a._control_matrix).sum()/cm.size, 0.9)
+
+
     def test_pulse_correlations(self):
         """Test calculating pulse correlation quantities."""
         for d, n_dt in zip(testutil.rng.integers(2, 7, 11),
diff --git a/tests/test_util.py b/tests/test_util.py
index f5003d6..98edbc3 100644
--- a/tests/test_util.py
+++ b/tests/test_util.py
@@ -559,7 +559,7 @@ def test_progressbar_range(self):
 
     def test_parse_optional_parameters(self):
 
-        @util.parse_optional_parameters({'foo': [1, 'bar', (2, 3)]})
+        @util.parse_optional_parameters(foo=[1, 'bar', (2, 3)], x=(2, 3))
         def foobar(a, b, foo=None, x=2):
             pass
 
@@ -581,6 +581,12 @@ def foobar(a, b, foo=None, x=2):
                              f"Invalid value for foo: {[1, 2]}." +
                              f" Should be one of {[1, 'bar', (2, 3)]}")
 
+        with self.assertRaises(ValueError):
+            foobar(1, 1, [1, 2], 4)
+            self.assertEqual(str(err.exception),
+                             f"Invalid value for x: {4}." +
+                             f" Should be one of {(2, 3)}")
+
 
 @pytest.mark.skipif(
     qutip is None,