Authored by M. Sifft and D. Hägele

We are excited to introduce MarkovAnalyzer, a high-speed Python toolkit designed for the analysis of hidden Markov models. This toolkit offers two distinct approaches: the well-established forward algorithm and our innovative polyspectra fitting method, both of which facilitate the study of complex Markov systems. The forward algorithm is a widely recognized technique, detailed on Wikipedia. In contrast, the polyspectra method is a novel contribution from our team at Ruhr University Bochum. This approach involves a comparative analysis between the experimental polyspectra—advanced extensions of the traditional power spectrum—and their theoretical models, which are derived from a given transition matrix and a measurement operator. This operator translates the state of the Markov system into measurable values. By matching the theoretical polyspectra with the experimental data, we can accurately estimate variable parameters within the hidden Markov model. The use of polyspectra in analyzing Markov processes comes with several significant benefits:

• Direct Input Utilization: Our method accepts the raw, experimentally measured Markov process data. This eliminates the need to categorize noisy measurements into distinct output levels, allowing for the analysis of data where noise is prevalent and distinct levels are obscured.
• Speed and Efficiency: Polyspectra-based analysis can outpace the forward algorithm, especially with large datasets. When dealing with terabyte-sized data, polyspectra need to be computed only once. These computations result in compact kilobyte-sized data, which are then used exclusively for the remainder of the analysis process.

MarkovAnalyzer is poised to revolutionize the way researchers approach the analysis of Markov systems, offering both efficiency and precision in one user-friendly package.

Installation

MarkovAnalyzer is available on pip and can be installed with

pip install markovanalyzer

Installation of Arrayfire

Besides running on CPU, MarkovAnalyzer offers GPU support for Nvidia and AMD cards. Depending on the hardware used, the usage of a GPU is highly recommended for Markov systems with more than about 100 states. A comprehensive installation guide for Linux + NVidia GPU can be found here. For GPU calculations the high performance library Arrayfire is used. The Python wrapper (see here) is automatically installed when installing SignalSnap, however, ArrayFire C/C++ libraries need to be installed separately. Instructions can be found can be found here and here.

Documentation

The documentation of the package can be found here. The package is divided into two parts: the polyspectra-calculator module, the fitting-tools module, and the forward module.

Polyspectra-Calculator Module

The Simulation Module allows for the calculation of the theoretical quantum polyspectra directly from the system’s transition matrix.

Fitting-Tools Module

The Fitting-Tools Module enables a user-friendly characterization of a Markov system in the lab based on the polyspectra of a measurement of a Markov process. These polyspectra can be calculated via our SignalSnap package. After providing a model transition matrix with one or more variable parameters, these parameters are estimated by fitting the theoretical model prediction of the polyspectra to their measured counterparts.

Forward Module

Here, the Forward Algorithm for the estimation of Markov models is implemented.

Example: Characterization of a Two-State Markov Model via Polyspectra

We want to deduce a Markov model (i.e., it’s transition matrix) from the observation of the Markov process. We are given data, that looks like this:

This is only a short excerpt of the full 6 min dataset. Using the SignalSnap library we are firstly calculating the polyspectra of that measurement. More details about the SignalSnap Code can be found on its GitHub page.

from signalsnap import SpectrumCalculator, SpectrumConfig, PlotConfig
import numpy as np
import h5py

Here is polyspectra are calculated and stored.

path = 'example_data/long_measurement.h5'
group_key = 'day1'
dataset = 'measurement1'

config = SpectrumConfig(dataset=dataset, group_key=group_key, path=path, f_unit='Hz', 
                        spectrum_size=150, f_max=2000, order_in=[1,2,3,4], 
                        backend='cpu')

spec = SpectrumCalculator(config)

f, s, serr = spec.calc_spec()

plot_config = PlotConfig(plot_orders=[2,3,4], arcsinh_plot=False, arcsinh_const=0.0002)
fig = spec.plot(plot_config)

path = 'example_data/two_state_example_spectra.pkl'
spec.save_spec(path, remove_S_stationarity=True)

Characterization is performed by fitting the theoretical polyspectra of a Markov model with variable parameters to their experimental counterparts calculated above. We are assuming a two-state model for the system that produced the data; hence, we begin with defining such a model with two variable transition rates.
The system undergoes transitions from the 0 to 1 state and from the 1 to the 0 state at rates gamma_01 and gamma_10, respectively. Each needs to be associated with a measurement value. From the histogram above we know that the value of the 0 state might be around 0, whereas the value of the 1 state might be around 26. Since we don’t know for sure, we leave these measurement values also as variable parameters n_0 and n_1.

import markovanalyzer as ma

def set_system(params):
      
    rates = {'0->1': params['gamma_01'],
             '1->0': params['gamma_10']}
             
    m_op = np.array([params['n_0'],params['n_1']])
    
    markov_system = ma.System(rates, m_op)

    return markov_system

Now, we can set start values and bounds for the fit of the parameters. A parameter c is always part of the fit and acconts for constant white noise outset in the power spectrum.

system_fit = ma.FitSystem(set_system)

parameter = {'gamma_01': [2.0396975e+04, 0, 1e5, True],
             'gamma_10': [1.0345057e+04, 0, 1e5, True],
             'n_0': [0, 0, 1e8, True],
             'n_1': [20, 0, 1e8, True],
             'c': [-2.7651282e-01 , -1e14, 1e14, True]}

path = 'example_data/two_state_example_spectra.pkl'

result = system_fit.complete_fit(path, parameter, 
                        method='least_squares', xtol=1e-6, ftol=1e-6, show_plot=True, fit_modus='order_based',
                        fit_orders=(1,2,3,4), beta_offset=False)

Authored by A. Ghorbanietemad and D. Hägele

PolyPyLive is designed to seamlessly integrate SignalSnap with streaming devices such as TimeTagger, thus enabling real-time data acquisition and processing. Moreover, PolyPyLive’s architecture is crafted with compatibility in mind, allowing for the direct analysis of collected data and saved data in SignalSnap without the need for any modifications. It ensures a smooth workflow transition from real-time processing to in-depth analysis, providing a comprehensive toolset for measurement.

Key Advantages

• Real-Time Analysis: PolyPyLive is unique in its ability to calculate and display higher-order spectra as data is being acquired.
• Advanced Estimation Techniques: The library employs SignalSnap, which utilizes an unbiased cumulant-based estimator for general signals and a moment-based estimator for coherent signals, ensuring accurate spectral computation.
• Enhanced Visibility: Users can opt for an arcsinh scale for signal display, improving the interpretability of results in real-time.

Testing and Readiness

PolyPyLive has been tested for detcting radio signals and is now ready for deployment in laboratory settings. As we strive to refine and enhance its functionality, we actively seek feedback from experimentalists. Your insights are invaluable to us for driving further improvements and ensuring that PolyPyLive meets the real-world demands of signal measurement research and development. We encourage users to share their experiences and suggestions. This revision emphasizes the readiness of the library for laboratory use while also making a clear call to action for feedback from the community, which is crucial for the iterative development process.

Caution and Future Development

This initial alpha version of PolyPyLive supports single-channel data acquisition. We are actively working to extend its capabilities to include multi-channel support in upcoming releases. Please be aware that the software saves data with the filename data.pkl by default; users have the flexibility to choose the storage directory path.

Example/TimeTagger Swabian Instrument

# This script can be used by cloning the repository and adjusting
# the parameters based on your needs.
# In the future, it will be developed into a Python package.
# ========================== Example usage ==========================
tagger_controller = _TaggerController()

channel_config = _ChannelConfig(channel=2,
                                trigger_level=1.5,
                                dead_time=0,
                                falling=False)

data_acq_config = _DataAcqConfig(buffer_size=10000000,
                                 time_unit='us',
                                 duration=int(30e7),
                                 delay=0,
                                 time_stamp_measure='s',
                                 save=False)

signal_config = _SignalConfig(signal_choice='S2',
                              backend='cuda',
                              m = 10,
                              f_max=4000,
                              coherent=False)

plot_config = _PlotConfig(data_points=200,
                          arcsinh_scale=False,
                          arcsinh_const=0.001,
                          sigma=3)

stream_setup = _StreamSetup(tagger_controller.tagger, channel_config, data_acq_config, signal_config, plot_config)

stream_setup.start()
stream_setup.process_data()  # Process and print the data
tagger_controller.free()

We can look at this code in details in the following:

tagger_controller = _TaggerController()

It is always needed since it communicates with the hardware and creates the TimeTagger object. If there is no signal generator, setTestSignal makes sure that you are still able to test the program by adding:

tagger_controller = _TaggerController()
tagger_controller.enable_test_signal([2])

to the program. The number 2 here is the number of channel in use.

Now it is time for the configuration of your data acquisition:

data_acq_config = _DataAcqConfig(buffer_size=10000000,
                                 time_unit='us',
                                 duration=int(30e7),
                                 delay=0,
                                 time_stamp_measure='s',
                                 save=True,
                                 save_path='your path')

Here ,you can set your buffer (higher buffer insures that your data are not being overwritten in TimeTagger). Then, you set the time unit for the duration of your measurement and delay as well as the duration. You also have to choose the time stamps measure. Saving path is your choice. In case of None, it will be saved on your desktop.

Next part of the program allows you to configurate your Channel.

channel_config = _ChannelConfig(channel=2,
                                trigger_level=1.5,
                                dead_time=0,
                                falling=False)

In the next step, you can verify your real-time visualization spectrum. This means that even if you select ‚S2‘, the power spectrum will be displayed. However, in the background, all ‚S1‘, ‚S2‘, ‚S3‘, and ‚S4‘, as well as their averages, are being calculated. You can select the backend according to signalsnap. The number mis the number of windows SignalSnaps needs to estimate cumulants. The larger m the faster the calculation. Furthermore, you need to specify a maximum frequency in Hz. For example, if you set f_max=2000, it means 2000Hz. If you choose a coherent signal, the estimators will be moment-based. If you select coherent=False, the estimators will be cumulant-based. It is optional to specify a minimum frequency. However, such setting automatically excludes ‚S3‘ from calculations.

signal_config = _SignalConfig(signal_choice='S2',
                              backend='cuda',
                              m = 10,
                              f_max=2000,
                              coherent=True)

In the end of the configuaraton, you make you plots to your liking:

plot_config = _PlotConfig(data_points=200,
                          arcsinh_scale=False,
                          arcsinh_const=0.001,
                          sigma=3)

The data_points parameter determines the resolution of your spectra. If you set arcsine_scale to True, it will ensure that small fluctuations are visible according to the arcsinh_const parameter. The smaller the value of this constant, the clearer the small fluctuations will be. The sigma parameter represents the significance value. It’s important to note that these configurations will NOT be saved in your saved data.

Additionally, the end part of the program is designed to allow the configurations to communicate with each other and to close the TimeTagger after the measurement is complete.

stream_setup = _StreamSetup(tagger_controller.tagger, channel_config, data_acq_config, signal_config, plot_config) stream_setup.start() stream_setup.process_data() tagger_controller.free()

Authored by M. Sifft and D. Hägele

The QuantumCatch package is open-source software for simulating and evaluating continuous quantum measurements via so-called quantum polyspectra. Here we refer to the polyspectra as second to fourth order spectra (powerspectrum, bispectrum, and 2D cut through trispectrum). The simulation of measurement traces (integration of the stochastic master equation) is implemented via the QuTiP toolbox whereas the calculation of polyspectra from Hamiltonians or measurements traces recorded in the lab is performed as described in this paper and this paper which also shows the utilization of quantum polyspectra to extract Hamiltonian parameters from a quantum measurement.

The quantum polyspectra approach enable the characterization of very general quantum systems that may include

• Environmental damping
• Measurement backaction (Zeno effect) and arbitrary measurement strength
• Coherent quantum dynamics
• Stochastic in- and out-tunneling
• Single-photon measurements
• Additional detector noise
• Simultaneous measurement of non-commuting observables
• Incorporation of temperatures
• Completely automatic analysis of arbitrary measurement traces
• Covers all limiting case of weak spin noise measurements, strong measurements resulting in quantum jumps, and single photon sampling

This poster provides an overview of the quantum polyspectra approach to quantum system characterization. Here is a brief summary: The analysis of a continuous measurement record 𝑧(𝑡) poses a fundamental challenge in quantum measurement theory. Different approaches have been used in the past as records can, e.g., exhibit predominantly Gaussian noise, telegraph noise, or clicks at random times. This poster summarizes our latest findings, showing that quantum measurements from all the aforementioned cases can be analyzed in terms of higher-order temporal correlations of the detector output 𝑧(𝑡) and be related to the Liouvillian of the measured quantum system. The comparison of temporal correlations via so-called quantum polyspectra is based on expressions derived without approximation from the stochastic master equation [1] and automated without requiring manual post-processing of the detector output. This allows for fitting of system parameters such as tunneling rates in a quantum transport experiment [2]. The very general stochastic master equation approach includes coherent quantum dynamics, environmental damping, and measurement backaction at arbitrary measurement strength. This enables a systematic evaluation of quantum measurements from the realms of conventional spin noise spectroscopy, quantum transport experiments, and ultra-weak measurements with stochastically arriving single photons [3,4]. [1] Hägele et al., PRB 98, 205143 (2018), [2] Sifft et al., PRR 3, 033123 (2021), [3] Sifft et al. PRA 107, 052203, [4] Sifft et al., arXiv:2310.10464

Installation

QuantumCatch is available on pip and can be installed with

pip install quantumcatch

Installation of Arrayfire

Besides running on CPU, QuantumCatch offers GPU support for Nvidia and AMD cards. Depending on the hardware used, the usage of a GPU is highly recommended for quantum systems with more than about 10 states. A comprehensive installation guide for Linux + NVidia GPU can be found here. For GPU calculations the high performance library Arrayfire is used. The Python wrapper (see here) is automatically installed when installing SignalSnap, however, ArrayFire C/C++ libraries need to be installed separately. Instructions can be found can be found here and here.

Documentation

The documentation of the package can be found here. The package is divided into two parts: the simulation module and the fitting-tools module.

Simulation Module

The Simulation Module allows for the integration of the Stochastic Master Equation via the QuTip library. This way an example measurement trace can be created. Moreover, it allows for the calculation of the theoretical quantum polyspectra directly from the system Liouvillian.

Fitting-Tools Module

The Fitting-Tools Module enables a user-friendly characterization of a quantum system in the lab based on the quantum polyspectra of an experimental continuous quantum measurement. These quantum polyspectra can be calculated via our SignalSnap package. After providing a model Liouvillian with one or more variable parameters, these parameters are estimated by fitting the theoretical model prediction of the polyspectra to their measured counterparts.

Example: Continuous Measurement of a Qunatum Dot’s Occupation

Here we are demonstrating how to simulate a simple quantum point contact measurement of the charge state of a quantum dot. First we have to import the QuantumCatch package. We will also import QuTip and NumPy. The analysis and generation module are imported automatically.

import quantumcatch as qc
import numpy as np
from qutip import *

Next, we are going to define the quantum dot system. The system has a Hamiltonian of H = 0 . The tunneling of electrons into and from the quantum dot and the measurement are modelled via Lindbladian damping terms while the measurement gets an additional backaction as defined in the stochastic master equation:

Let’s see how this system is implemented for QuantumCatch using QuTip functions. First we are going to define the operators:

dot_levels = 2
a = destroy(dot_levels)
n = a.dag() * a
rho_0 = ket2dm(basis(dot_levels, 0)) + ket2dm(basis(dot_levels, 1)) 
rho_0 /= np.trace(rho_0)
H = 0 * n

Then we are going to define the initial density matrix of the quantum dot. Here we are starting in a superposition of full and empty:

rho_0 = ket2dm(basis(dot_levels, 0)) + ket2dm(basis(dot_levels, 1)) 
rho_0 /= np.trace(rho_0)

By defining the tunneling rate and measurement strength we are also setting the units of time. Here we are choosing ms so the rates will be given in kHz.

beta = 5  # in sqrt(kHz)
gamma_out = 0.5  # in kHz
gamma_in = 0.5  # in kHz

Now we have to add damping (c_ops) and measurement terms (sc_ops). We begin by defining the damping/measurement strength γ and β choosing the name of the corresponding operator as key an the dictionary:

c_measure_strength = {
    'a': gamma_out**0.5,
    'a_d': gamma_in**0.5,
}

sc_measure_strength = {
    'n': beta
}

Afterward we define the damping and measurement operator in QuTiP fashion:

c_ops = {
    'a': a,
    'a_d': a.dag(),
}

sc_ops =  {
    'n': n,
}

For the calculation of the observable, we are also defining the operators with which the expectation values are calculated:

e_ops = {
    'n': n,
}

The System is now ready to be initialized as follows:

system = qc.System(H, rho_0, c_ops, sc_ops, e_ops, c_measure_strength, sc_measure_strength)

To run the simulation we simply have to define the desired time steps. Note that especially during a strong measurement narrow time steps should be chosen. The calc_trnasient function has many options. Please refer to the docs.

t_start = 0
t_stop = int(4e1)
n_steps = int(10e3)
t = np.linspace(t_start, t_stop, n_steps)
res = system.calc_transient(t, seed=None, _solver='milstein', _normalize=False)

The result can be visualized by providing a list of view operators as follows:

view_ops = {'n':1}
system.plot_transient(view_ops)

The corresponding quantum polyspectra can be calculated easily via calculate_spectrum function. Refer to the documentation for the option of this function. For example, here the usage of the GPU can be enabled setting enable_gpu=True.

f_start = 0
f_stop = 2
n_points = 200
fs = np.linspace(f_start, f_stop, n_points)

measure_op = 'n'

spec = system.calculate_spectrum(fs, order_in=[1,2,3,4], measure_op=measure_op)

The polyspectra can be plotted with the plot method:

fig = system.plot()

Example: From Quantum Polyspectra to the Model Liouvillian

We begin with formulating a model quantum system using standard QuTiP functions.

import quantumcatch as qc
import numpy as np
from qutip import *

The system has to by defined as a function (here set_system) that takes in an array of variable parameters. We will define the parameters soon, here you simply use the ´params´ dictionary with you own keys of choice (here gamma_in, gamma_out, and beta). These parameters will be determined by the fitting routine.

def set_system(params):
    # ------ Operators and Hamiltonian -----
    dot_levels = 2
    a = destroy(dot_levels)
    n = a.dag() * a
    rho_0 = ket2dm(basis(dot_levels, 0)) + ket2dm(basis(dot_levels, 1)) 
    rho_0 /= np.trace(rho_0)
    H = 0 * n
    
    # ------ System Parameters ------
    c_measure_strength = {
        'a': params['gamma_out']**0.5,
        'a_d': params['gamma_in']**0.5,
    }
    
    sc_measure_strength = {
        'n': params['beta']
    }
    
    c_ops = {
        'a': a,
        'a_d': a.dag(),
    }
    
    sc_ops =  {
        'n': n,
    }
    
    e_ops = {
        'n': n,
    }

    system = qc.System(H, rho_0, c_ops, sc_ops, e_ops, c_measure_strength, sc_measure_strength)
    
    return system, sc_ops, sc_measure_strength

With this parameterize system and the name of our measurement operator, we now create a fit object:

m_op = 'n'
system_fit = qc.FitSystem(set_system, m_op)

Now, we can define the parameters with their starting values and bounds in lmfit fashion. Assuming that the path leads to spectra calculated and saved with the SignalSnap library, the fit routine is started with the complete_fit function. Refer to the documentation for all possible function parameters.

parameter = {'gamma_in':[300, 0, 1e5, True],
             'gamma_out': [700, 0, 1e5, True],
             'beta': [1, 0, 1e12, True],
             'c': [0, -1e12, 1e12, True]}
       
path = 'example_data/polyspectra_for_fit_example.pkl'         
out = system_fit.complete_fit(path, parameter, method='least_squares', xtol=1e-8, ftol=1e-8, fit_modus='resolution_based',
                                                            fit_orders=(1,2,3,4))

This figure presents a comprehensive set of nine plots designed to evaluate the quality of the fit. The plots are arranged in a 3×3 grid and are interpreted as follows:

• (0,0): Depicts the first-order spectrum, S1, which is proportional to the mean of the measurement. The experimental S1 trace is displayed in blue, accompanied by its 3σ (three standard deviations) error. The fitted S1 is overlaid in orange.
• (0,1): Illustrates both the experimental and theoretical power spectra.
• (0,2): Presents the relative error between the experimental and theoretical power spectra, normalized to the experimental power spectrum, as a blue line. This line can be compared to the 3σ error of the power spectra values.
• (1,0): Displays the theoretical (upper left triangle) and experimental (lower right triangle) bispectra. This combined representation is facilitated by the symmetry properties of the bispectrum, with the diagonal serving as a symmetry axis.
• (1,1): Shows the relative error between the experimental and theoretical bispectra, normalized to the experimental bispectrum. A pixel is colored green when the difference between the experimental and theoretical bispectra is less than the 3σ error of the spectral values.
• (1,2): Features experimental and theoretical bispectral values along the axes (ω1, 0) and the diagonal, accompanied by a 3σ error band.
• (2,0) to (2,2): These plots mirror the type of plots shown for the bispectrum, but now with the trispectrum.

Each plot provides a unique perspective on the data, collectively offering a robust assessment of the fit’s accuracy.

Support

The development of the QunatumCatch package is supported by the working group Spectroscopy of Condensed Matter of the Faculty of Physics and Astronomy at the Ruhr University Bochum.

Dependencies

For the package multiple libraries are used for the numerics and displaying the results:

IPyWidgets

Numpy

Matplotlib

H5py

ArrayFire

Numba

SciPy

Tqdm

Plotly

Psutil

Cachetools

Signalsnap

Pandas

Lmfit

Authored by M. Sifft, A. Ghorbanietemad and D. Hägele

We present a fast Python toolbox for higher-order spectral analysis of time series. The usual second-order power spectrum and its higher-order generalization – so-called bi- and trispectra – are efficiently calculated on any platform. The toolbox supports GPU acceleration using the ArrayFire library. The treatment of large datasets is not limited by available RAM. We were able to process 11.6 GB of data (given in the hdf5 format) within just one minute to obtain a power spectrum with a one-million-point resolution using a five-year-old Nvidia GeForce GTX 1080. Similarly, 1000×1000 points of a trispectrum were obtained from processing 3.3 GB of data per minute.

Here are a few outstanding features of SignalSnap:

• Errors of spectral values are automatically calculated (beginners example, example)
• Support for just-in-time import from hdf data (dataset does not have to fit in RAM) (example)
• Functions for conversion of Numpy array to hdf data is also provided (example)
• Functions for storing and loading calculated spectra together with metadata (example)
• Correlations between two time series can be calculated (example)
• All calculation can be performed on GPU (NVidia and AMD) (see Arrayfire) (example)
• Advanced plotting options for two-dimensional higher-order spectra (as seen in most examples)
• Usage of unbiased estimators for higher-order cumulants (see Literature below)
• Efficient implementation of the confined Gaussian window for an optimal RMS time-bandwidth product (see Literature below)
• Special functions for the verification of the stationarity of a signal (example)
• Spectra can be calculated from timestamps instead of a continuous trace (example)

Installation

SignalSnap is available on pip and can be installed with

pip install signalsnap

Installation of Arrayfire

A comprehensive installation guide for Linux + NVidia GPU can be found here. For GPU calculations the high performance library Arrayfire is used. The Python wrapper (see here) is automatically installed when installing SignalSnap, however, ArrayFire C/C++ libraries need to be installed separately. Instructioins can be found can be found here and here.

Documentation

A documentation of SignalSnap’s functions can be found here.

Examples

Examples for every function of the package are currently added to the folder Examples. Here are a few lines to get you started. We will generate some white noise as signal/dataset and store it as Numpy array called y.

from signalsnap import SpectrumCalculator, SpectrumConfig, PlotConfig
import numpy as np

rng = np.random.default_rng()

# ------ Generation of white noise -----
f_unit = 'kHz'
fs = 10e3  # sampling rate in kHz
N = 1e5  # number of points
t = np.arange(N) / fs  # unit is automatically chosen to be 1/f_unit = ms
y = rng.normal(scale=1, size=t.shape)

Now we creat a spectrum object and feed it with the data. This object will store the dataset, later the spectra and errors, all freely chosen variables and contains the methods for calculating the spectra, plotting and storing.

config = SpectrumConfig(data=y, delta_t=1/fs, f_unit='kHz', 
                        spectrum_size=128, order_in=[2], 
                        f_max=5e3, backend='cpu')
spec = SpectrumCalculator(config)
f, s, serr = spec.calc_spec()

T_window: 2.540e-02 ms
Maximum frequency: 5.000e+03 kHz

The output will show you the actual length of a window (in case your T_window is not a multiple of 1/fs), the maximum frequency (Nyquist frequency) and the number of point of the calculated spectrum. The data points in the first window are plotted, so you can verify the window length (which is also given in points by chunk shape). The function will return f the frequencies at which the spectrum has been calculated, s the spectral values, and serr the error of the spectra value (1 sigma).

Visualization of the results is just as easy as the calculation.

plot_config = PlotConfig(plot_orders=[2], plot_f_max=5e3/2)
fig = spec.plot(plot_config)

Besides the power spectrum (blue line) the error bands (1 to 5 sigma) are shown as grey lines in the plot. Now, we can even verify that we are dealing with true Gaussian noise by calculating the higher-order spectra of the time series.

Why higher-order spectra?

Higher-order spectra contain additional information that is not contained within a power spectrum. The toolbox is capable of calculating the third- and four-order spectrum (also called bi- and trispectrum, respectively). These have the following properties:

• Spectra beyond second order are not sensitive to Gaussian noise.
• Bispectrum: shows contributions whenever the phase of two frequencies and their sum are phase correlated (e.g. by mixing two signals)
• Trispectrum: can be interpreted as intensity correlation between two frequencies

Let’s calculate all spectra up to fourth order of the dataset above and verify that the signal does not deviate significantly from Gaussian noise using the first property (has no significant higher-order contributions). We only have to change the order_in argument:

config = SpectrumConfig(data=y, delta_t=1/fs, f_unit='kHz', 
                        spectrum_size=128, order_in='all', 
                        f_max=5e3, backend='cpu')
spec = SpectrumCalculator(config)
f, s, serr = spec.calc_spec()

Plotting can also be done as before by changing the order_in argument:

plot_config = PlotConfig(plot_orders=[2,3,4], plot_f_max=5e3/2, green_alpha=0)
fig = spec.plot(plot_config)

Now, the third-and fourth-order spectra (S3 and S4) are visible. Just like the power spectrum they are noisy. To decide which of the fluctuations are significant we need a way of displaying errors in the two-dimensional plots. Here, errors are visualized be overlaying a green color on the spectral contributions which deviate from zero less than a certain number of standard deviations.

plot_config = PlotConfig(plot_orders=[2,3,4], plot_f_max=5e3/2, sigma=3)
fig = spec.plot(plot_config)

Clearly, all higher-order contributions are nothing but noise and we have, therefore, verifed that our original dataset was Gaussian noise (and even white noise due to the flat power spectrum).

Support

The development of the SignalSnap package is supported by the working group Spectroscopy of Condensed Matter of the Faculty of Physics and Astronomy at the Ruhr University Bochum.

Dependencies

For the package multiple libraries are used for the numerics and displaying the results:

• NumPy
• SciPy
• MatPlotLib
• tqdm
• Numba
• h5py
• ArrayFire

Literature

Unbiased estimators are used for the calculation of higher-order cumulants. Their derivation can be found in this paper. An explanation for how the spectra are calculated can be found in Appendix B of this paper.

Authored by S. Starosielec and D. Hagele

We provide easy to use MATLAB code (CGWindow.zip) for calculating the Confined Gaussian Window (CGW) and Approximate Confined Gaussian Window (ACGW) which posses an optimal RMS time-bandwidth product. A single command w = replacewinbyCG(w); replaces a suboptimal window w from existing code by the Confined Gaussian window. The new window conserves the RMS temporal width and the square norm w’*w of the initial window.