diff --git a/src/python_spectrometer/__init__.py b/src/python_spectrometer/__init__.py
index a6319b875357cefede1bf2eecefedc201e98eadc..5b54107ecac6dcd05110c952089def4a3ec3bf17 100644
--- a/src/python_spectrometer/__init__.py
+++ b/src/python_spectrometer/__init__.py
@@ -93,8 +93,11 @@ Finally, plot options can be changed dynamically at runtime::
spect.plot_raw = True # Updates the figure accordingly
spect.plot_timetrace = False
+Examples
+--------
+
Example from :func:`scipy.signal.welch`
----------------------------------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
In this short demonstration, we reproduce the example from
:func:`scipy:scipy.signal.welch`. To this end, we write a custom
``DAQ`` class that generates a noisy sine signal.
@@ -146,6 +149,46 @@ Finally, we can also plot data in dB relative to a given dataset.
>>> float(data.max()) # Factor two in amplitude is approx 3 dB
3.0284739712568682
+Analyzing filter behavior
+^^^^^^^^^^^^^^^^^^^^^^^^^
+:mod:`qutil:qutil.signal_processing.real_space` and
+:mod:`qutil:qutil.signal_processing.fourier_space` define filters that
+work in the time- and frequency-domain, respectively. We can visualize
+the filter properties using the spectrometer:
+
+>>> from tempfile import mkdtemp
+>>> import qutil.signal_processing as sp
+>>> from qutil.functools import partial
+>>> from python_spectrometer import daq, Spectrometer
+
+>>> def compare_filters(typ: str, order: int):
+... spect = Spectrometer(daq.QoptColoredNoise(), savepath=mkdtemp(),
+... plot_dB_scale=True, plot_density=False,
+... threaded_acquisition=False)
+... spect.take('Baseline', n_seg=10, fs=1e4, df=0.1)
+... spect.procfn = getattr(sp.real_space, f'{typ}_filter')
+... spect.take(f'Real space {order}. order {typ} filter',
+... n_seg=10, f_max=1e2, fs=1e4, df=0.1, order=order)
+... spect.procfn = sp.real_space.Id
+... spect.psd_estimator = partial(
+... sp.real_space.welch,
+... fourier_procfn=getattr(sp.fourier_space, f'{typ}_filter')
+... )
+... spect.take(f'Fourier space {order}. order {typ} filter',
+... n_seg=10, f_max=1e2, fs=1e4, df=0.1, order=order)
+
+RC and Butterworth first order filters are the same (up to real-space
+implementation):
+
+>>> compare_filters('RC', 1)
+>>> compare_filters('butter', 1)
+
+For higher orders, they differ:
+
+>>> compare_filters('RC', 5)
+>>> compare_filters('butter', 5)
+
+
See the documentation of :class:`~core.Spectrometer` and its methods
for more information.
"""
diff --git a/src/python_spectrometer/_plot_manager.py b/src/python_spectrometer/_plot_manager.py
index 24030f2eb792fc8b7b98f76f676a3f571a48ae6d..97865a8f284a567c557231df697b2d696ef7fe7f 100644
--- a/src/python_spectrometer/_plot_manager.py
+++ b/src/python_spectrometer/_plot_manager.py
@@ -74,10 +74,6 @@ class PlotManager:
self.legend_kw.setdefault('loc', 'upper right')
self.figure_kw.setdefault('num', f'Spectrometer {len(self.__instances) + 1}')
- if not any('layout' in key for key in self.figure_kw.keys()):
- # depending on the matplotlib version, this setting is either
- # layout='tight' or tight_layout=True
- self.figure_kw.setdefault('layout', 'tight')
if self.subplot_kw.pop('sharex', None) is not None:
warnings.warn('sharex in subplot_kw not negotiable, dropping', UserWarning)
@@ -102,7 +98,7 @@ class PlotManager:
try:
self._fig = plt.figure(**self.figure_kw)
except TypeError:
- if layout := self.figure_kw.pop('layout', False):
+ if layout := self.figure_kw.pop('layout', None) is not None:
# matplotlib < 3.5 doesn't support layout kwarg yet
self.figure_kw[f'{layout}_layout'] = True
elif layout is False:
@@ -579,6 +575,13 @@ class PlotManager:
if self._leg is not None:
self._leg.remove()
+ if 'layout' not in self.figure_kw:
+ try:
+ self.fig.set_layout_engine('tight')
+ self.fig.get_layout_engine().execute(self.fig)
+ finally:
+ self.fig.set_layout_engine('none')
+
self.fig.canvas.draw_idle()
self.fig.canvas.flush_events()
diff --git a/src/python_spectrometer/core.py b/src/python_spectrometer/core.py
index ec130ba01dfdf244528339bc277db7bfecc0838d..b0d4703f6c10a23bd4cab9c660e5ecbc6296fffc 100644
--- a/src/python_spectrometer/core.py
+++ b/src/python_spectrometer/core.py
@@ -1089,6 +1089,8 @@ class Spectrometer:
# TODO: Instead of setting up separate LiveView's for time and frequency
# data, the cleaner option would arguably be to subclass and add
# another subplot.
+ # TODO: When plot_timetrace=True and one view is paused, closing a
+ # figure causes the thread to not join.
self._assert_ready()
@@ -1109,7 +1111,10 @@ class Spectrometer:
t = np.linspace(0, T, settings['n_pts'])
while not stop_event.is_set():
*_, yt = get_queue.get()
- put_queue.put((t, yt))
+ if np.issubdtype(self.daq.DTYPE, np.complexfloating):
+ put_queue.put((t, {'X': yt.real, 'Y': yt.imag}))
+ else:
+ put_queue.put((t, yt))
get_queue.task_done()
live_view.flush_queue(get_queue)
@@ -1153,7 +1158,7 @@ class Spectrometer:
settings = self.daq.setup(**settings)
T = settings['n_pts'] / settings['fs']
- if np.issubdtype(self.daq.DTYPE, np.complexfloating):
+ if np.issubdtype(self.daq.DTYPE, np.complexfloating) and self.plot_negative_frequencies:
freq_xscale = _asinh_scale_maybe()
xlim = np.array([-settings['f_max'], settings['f_max']])
else:
@@ -1167,7 +1172,11 @@ class Spectrometer:
live_view_kw.setdefault('blocking_queue', True)
live_view_kw.setdefault('autoscale', 'c')
- live_view_kw.setdefault('autoscale_interval', None)
+ live_view_kw.setdefault('autoscale_interval_ms', None)
+ live_view_kw.setdefault(
+ 'plot_legend',
+ 'upper right' if np.issubdtype(self.daq.DTYPE, np.complexfloating) else False
+ )
fixed_kw = dict(
plot_line=True, xlim=xlim, ylim=(0, (max_rows - 1) * T),
@@ -1187,6 +1196,7 @@ class Spectrometer:
if self.plot_timetrace:
fixed_kw = dict(xlim=(0, T), xlabel='$t$', ylabel='Amplitude',
+ n_lines=2 if np.issubdtype(self.daq.DTYPE, np.complexfloating) else 1,
units={'x': 's', 'y': self.processed_unit})
time_kw = live_view_kw | fixed_kw
if any(time_kw[key] != val for key, val in live_view_kw.items()):
diff --git a/src/python_spectrometer/daq/__init__.pyi b/src/python_spectrometer/daq/__init__.pyi
index 1fd2d38ff8846ee49d3ad65731b147d21335b995..83c91696a8ddb3effe78cc899bb549e0ec6aad4a 100644
--- a/src/python_spectrometer/daq/__init__.pyi
+++ b/src/python_spectrometer/daq/__init__.pyi
@@ -1,5 +1,5 @@
-__all__ = ['atsaverage', 'atssimple', 'qcodes', 'simulator', 'swabian_instruments', 'zurich_instruments',
- 'AlazarATS9xx0', 'Keysight344xxA', 'NationalInstrumentsUSB6003',
+__all__ = ['atsaverage', 'atssimple', 'qcodes', 'simulator', 'swabian_instruments',
+ 'zurich_instruments', 'AlazarATS9xx0', 'Keysight344xxA', 'NationalInstrumentsUSB6003',
'SwabianInstrumentsTimeTagger', 'DAQSettings', 'QoptColoredNoise',
'ZurichInstrumentsMFLIScope', 'ZurichInstrumentsMFLIDAQ']
diff --git a/src/python_spectrometer/daq/simulator.py b/src/python_spectrometer/daq/simulator.py
index ed8ebc179a4c6a793133a422a629b6e9506344f8..8e7915a88d750b393a5e97ab19d423f1aa22b442 100644
--- a/src/python_spectrometer/daq/simulator.py
+++ b/src/python_spectrometer/daq/simulator.py
@@ -1,29 +1,40 @@
"""Provides a simulation backend for data acquisition.
+The :meth:`DAQ.acquire` method of DAQs in this module accepts a
+``delay`` keyword argument which introduces a delay in between yields
+to simulate a finite data acquisition time. If True, delays by the
+amount of time it would take to actually acquire data with the given
+settings; if float delays by the given amount.
+
Examples
--------
>>> import python_spectrometer as pyspeck
>>> import tempfile
>>> speck = pyspeck.Spectrometer(pyspeck.daq.QoptColoredNoise(),
... savepath=tempfile.mkdtemp())
->>> speck.take('a test', fs=10e3) #doctest: +ELLIPSIS
+>>> speck.take('a test', fs=10e3)
...
>>> speck.block_until_ready() # for doctest
Add an artificial time delay to mimick finite data acquisition time:
->>> speck.take('delayed', n_avg=3, delay=True) #doctest: +ELLIPSIS
+>>> speck.take('delayed', n_avg=3, delay=True)
...
"""
from __future__ import annotations
+import copy
import dataclasses
+import inspect
import sys
import time
-from typing import Callable
+from collections.abc import Callable
+from typing import Literal, Union
import numpy as np
-from qutil.functools import partial
+from qutil.functools import partial, wraps
+from qutil.math import cexp
+from qutil.signal_processing import real_space
from .base import DAQ, AcquisitionGenerator
@@ -43,6 +54,64 @@ except ImportError as e:
"'pip install qopt.'") from e
+def with_delay(meth):
+ """Wraps an acquisition generator to accept the *delay* kwarg."""
+
+ @wraps(meth)
+ def wrapped(self, *, delay=False, **settings):
+ skip_delay = settings.pop('_skip_delay', False)
+
+ if delay is True:
+ delay = settings['n_pts'] / settings['fs']
+
+ it = meth(self, _skip_delay=True, **settings)
+ while True:
+ tic = time.perf_counter()
+ try:
+ data = next(it)
+ except StopIteration as stop:
+ return stop.value
+ else:
+ if delay and not skip_delay:
+ time.sleep(max(0, delay - (time.perf_counter() - tic)))
+
+ yield data
+
+ # Insert parameter sig
+ delay_param = inspect.Parameter('delay', inspect.Parameter.KEYWORD_ONLY, default=True,
+ annotation=Union[bool, float])
+
+ parameters = list(inspect.signature(meth).parameters.values())
+ if parameters[-1].kind is inspect.Parameter.VAR_KEYWORD:
+ parameters = parameters[:-1] + [delay_param, parameters[-1]]
+ else:
+ parameters = parameters + [delay_param]
+
+ wrapped.__signature__ = inspect.signature(wrapped).replace(parameters=parameters)
+ return wrapped
+
+
+class MonochromaticNoise(DAQ):
+ """Generate monochromatic sinusoidal noise with random phase.
+
+ This DAQ implementation produces sinusoidal data with a fixed frequency
+ but random phase for each acquisition, simulating a simple signal with noise.
+
+ Inherits from the base DAQ class and implements the required acquire method.
+ """
+
+ @with_delay
+ def acquire(self, *, n_avg: int, fs: float, n_pts: int, A: float = 1, f_0: float = 50,
+ **settings) -> AcquisitionGenerator[DAQ.DTYPE]:
+ """Generate sinusoidal data with random phase."""
+
+ t = np.arange(0, n_pts / fs, 1 / fs)
+ rng = np.random.default_rng()
+
+ for _ in range(n_avg):
+ yield np.sin(2 * np.pi * (t * f_0 + rng.random()))
+
+
@dataclasses.dataclass
class QoptColoredNoise(DAQ):
"""Simulates noise using :mod:`qopt:qopt`.
@@ -52,11 +121,17 @@ class QoptColoredNoise(DAQ):
:class:`~python_spectrometer.daq.settings.DAQSettings`
for more information on setup parameters.
+ Attributes
+ ----------
+ spectral_density : Callable[[NDArray, ...], NDArray]
+ A function that generates the power spectral density for given
+ frequencies. Defaults to white noise with scale parameter
+ ``S_0``.
+
See Also
--------
:func:`qopt:qopt.noise.fast_colored_noise`
For information on the simulation.
-
"""
spectral_density: Callable[[NDArray, ...], NDArray] = dataclasses.field(
default_factory=lambda: QoptColoredNoise.white_noise
@@ -74,35 +149,120 @@ class QoptColoredNoise(DAQ):
"""White noise power spectral density with amplitude S_0."""
return np.full_like(f, S_0)
- def acquire(self, *, n_avg: int, fs: float, n_pts: int, delay: bool | float = False,
+ @with_delay
+ def acquire(self, *, n_avg: int, fs: float, n_pts: int,
**settings) -> AcquisitionGenerator[DAQ.DTYPE]:
- """Executes a measurement and yields the resulting timetrace.
-
- Optionally set a delay to simulate a finite data acquisition
- time. If True, delays by the amount of time it would take to
- actually acquire data with the given settings; if float delay
- by the given amount.
- """
- if delay is True:
- delay = 1 / settings['df']
-
+ """Executes a measurement and yields the resulting timetrace."""
for _ in range(n_avg):
- tic = time.perf_counter()
- data = qopt.noise.fast_colored_noise(
+ yield qopt.noise.fast_colored_noise(
partial(
settings.get('spectral_density', self.spectral_density),
**settings
),
dt=1/fs, n_samples=n_pts, output_shape=()
)
- if delay:
- time.sleep(delay - (time.perf_counter() - tic))
- yield data
# This is the place to return metadata (possibly obtained from the instrument)
return {'qopt version': qopt.__version__}
+class DemodulatorQoptColoredNoise(QoptColoredNoise):
+ """Simulates demodulated noisy data for lock-in measurements.
+
+ Extends QoptColoredNoise to demodulate the simulated signal using
+ complex IQ-demodulation, similar to a lock-in amplifier. This
+ provides a realistic simulation of demodulated signals as would be
+ measured in experiments using lock-in amplification techniques.
+ """
+ DTYPE = np.complexfloating
+
+ @staticmethod
+ def demodulate(signal: np.ndarray, IQ: np.ndarray, **settings) -> np.ndarray:
+ """Demodulate signal using the provided IQ reference.
+
+ Performs complex demodulation by multiplying the signal with
+ the IQ reference and applying an RC filter.
+
+ Parameters
+ ----------
+ signal :
+ Input signal to demodulate
+ IQ :
+ Complex IQ reference for demodulation
+ **settings :
+ Settings for RC filter, including filter parameters
+
+ """
+ # Don't highpass filter
+ settings = copy.deepcopy(settings)
+ settings.pop('f_min', None)
+ return real_space.RC_filter(signal * IQ, **settings)
+
+ @with_delay
+ def acquire(self, *, n_avg: int, freq: float = 50, filter_order: int = 1,
+ modulate_signal: bool = False,
+ filter_method: Literal['forward', 'forward-backward'] = 'forward',
+ **settings) -> AcquisitionGenerator[DTYPE]:
+ r"""Simulate demodulated noisy data.
+
+ Generates simulated data and performs IQ demodulation, mimicking
+ the behavior of a lock-in amplifier. Can simulate either just
+ input noise or noise in the full signal path.
+
+ See [1]_ for an introduction to Lock-in amplification.
+
+ Parameters
+ ----------
+ n_avg :
+ Number of outer averages.
+ freq :
+ Modulation frequency.
+ filter_order :
+ RC filter order used to filter the demodulated signal.
+ modulate_signal :
+ Add the simulated noise to the modulation signal to mimic
+ noise picked up by a Lock-In signal travelling through some
+ DUT. Otherwise, mimics the noise at the input of the
+ amplifier.
+
+ In other words, simulate a Lock-In output connected to an
+ input, or just simulate the input.
+
+ Note that if True, noise is assumed to be additive, that
+ is,
+
+ .. math::
+
+ x(t) = s(t) + \delta(t)
+
+ with :math:`s(t)` the output signal and :math:`\delta(t)`
+ the noise.
+ filter_method :
+ See :func:`~qutil:qutil.signal_processing.real_space.RC_filter`.
+
+ Yields
+ ------
+ data :
+ Demodulated data in complex IQ-representation.
+
+ References
+ ----------
+ .. [1] https://www.zhinst.com/europe/en/resources/principles-of-lock-in-detection
+
+ """
+ t = np.arange(0, settings['n_pts'] / settings['fs'], 1 / settings['fs'])
+ # demodulation by √2 exp(-iωt) (ZI convention)
+ IQ = np.sqrt(2) * cexp(-2 * np.pi * freq * t)
+
+ yield from (
+ self.demodulate(IQ.real + data if modulate_signal else data, IQ,
+ order=filter_order, method=filter_method, **settings)
+ for data in super().acquire(n_avg=n_avg, **settings)
+ )
+
+ return {'qopt version': qopt.__version__}
+
+
@deprecated("Use QoptColoredNoise instead")
class qopt_colored_noise(QoptColoredNoise):
...
diff --git a/src/python_spectrometer/daq/swabian_instruments.py b/src/python_spectrometer/daq/swabian_instruments.py
index b6adfa77dd323eaf0d1b137664563fbb0e82e437..5ff94f76a828f473a0160f9e430b01f6eb1d4f11 100644
--- a/src/python_spectrometer/daq/swabian_instruments.py
+++ b/src/python_spectrometer/daq/swabian_instruments.py
@@ -53,7 +53,7 @@ class SwabianInstrumentsTimeTagger(DAQ):
_ = atexit.register(TimeTagger.freeTimeTagger, tagger)
spect = Spectrometer(daq.swabian_instruments.SwabianInstrumentsTimeTagger(tagger, [1, 2]),
- raw_unit='cts', processed_unit='cts')
+ raw_unit='cts')
"""
tagger: tt.TimeTagger = dataclasses.field()
diff --git a/src/python_spectrometer/daq/zurich_instruments.py b/src/python_spectrometer/daq/zurich_instruments.py
index 876fd1f64d95327030e7464da5f043f87eac581f..64d900573c81900095b4737d58a41cd8fdbe0b04 100644
--- a/src/python_spectrometer/daq/zurich_instruments.py
+++ b/src/python_spectrometer/daq/zurich_instruments.py
@@ -39,7 +39,15 @@ Compare their results::
spect_scope.take(n_pts=2**14, fs=14.6e3)
The DAQ spectrum should show a peak at :math:`f=-500\,\mathrm{Hz}`,
-corresponding to the oscillator frequency.
+corresponding to the oscillator frequency. This is the shifted 0 Hz
+peak of the instrument's $1/f$ noise.
+
+Use ``procfn`` to compute the phase noise spectrum using the DAQ module::
+
+ spect_daq.procfn = lambda x, **_: np.angle(x)
+ spect_daq.processed_unit = 'rad'
+ spect_daq.drop('all')
+ spect_daq.take()
"""
from __future__ import annotations
@@ -57,6 +65,8 @@ from packaging import version
from qutil.domains import DiscreteInterval, ExponentialDiscreteInterval
from scipy.special import gamma
from zhinst import toolkit
+from zhinst.core.errors import SampleLossError
+from zhinst.toolkit.exceptions import ToolkitError
from .base import DAQ, AcquisitionGenerator
from .settings import DAQSettings
@@ -78,9 +88,6 @@ if version.parse(toolkit.__version__) < version.parse('0.5.0'):
raise ImportError('This DAQ requires zhinst-toolkit >= 0.5.0. '
"You can install it by running 'pip install zhinst-toolkit>=0.5.0'.")
-from zhinst.core.errors import SampleLossError
-from zhinst.toolkit.exceptions import ToolkitError
-
logger = logging.getLogger(__name__)
ZhinstDeviceT = TypeVar('ZhinstDeviceT', bound=toolkit.driver.devices.base.BaseInstrument)
@@ -116,14 +123,6 @@ class ZurichInstrumentsMFLIDAQ(_ZurichInstrumentsDevice):
components, X + iY, and therefore the resulting spectrum is two-
sided.
- .. note::
-
- For measurements where the input is not connected to the output
- of the Lock-in, this means there is a coherent signal at the
- oscillator frequency present in the time series data. Due to
- non-idealities, its peak in the spectrum has a finite width and
- may thus overshadow spectral features present in the input.
-
Parameters
----------
session : toolkit.session.Session
@@ -171,28 +170,14 @@ class ZurichInstrumentsMFLIDAQ(_ZurichInstrumentsDevice):
return MFLIDAQSettings
- def setup(self, bandwidth: Union[str, float] = 'auto', filter_order: Optional[int] = None,
- freq: float = 50, **settings: Mapping) -> Dict[str, Any]:
+ def setup(self, filter_order: Optional[int] = None, freq: float = 50,
+ **settings: Mapping) -> Dict[str, Any]:
r"""Sets up the daq module to acquire time series data.
See [1]_ for information on lock-in measurements.
Parameters
----------
- bandwidth : Union[str, float], optional
- The demodulator noise-equivalent power (NEP) bandwidth.
- The default is 'auto', in which case it is set to f_max/2.
-
- The bandwidth is related to the time constant of the RC
- filter by
-
- .. math::
-
- \tau = \frac{\Gamma\left(n - \frac{1}{2}\right)}
- {4\sqrt{\pi}f_\mathrm{NEP}\Gamma(n)},
-
- where :math:`n` is the filter order
-
filter_order : int, optional
The filter order. Not changed if not given.
freq : float, optional
@@ -211,9 +196,23 @@ class ZurichInstrumentsMFLIDAQ(_ZurichInstrumentsDevice):
**settings : Mapping
Additional settings for data acqusition.
+ Notes
+ -----
+ The demodulator 3 dB bandwidth is chosen as $f_\text{max}$. The
+ noise-equivalent power (NEP) bandwidth is related to the time
+ constant of the RC filter by
+
+ .. math::
+
+ \tau = \frac{\Gamma\left(n - \frac{1}{2}\right)}
+ {4\sqrt{\pi}f_\mathrm{NEP}\Gamma(n)},
+
+ where :math:`n` is the filter order. By default, it is set to
+ ``fs/4``.
+
Raises
------
- RuntimeError
+ RuntimeError :
If settings are incompatible with the hardware.
Returns
@@ -226,18 +225,22 @@ class ZurichInstrumentsMFLIDAQ(_ZurichInstrumentsDevice):
.. [1] https://www.zhinst.com/europe/en/resources/principles-of-lock-in-detection
"""
- settings = self.DAQSettings(freq=freq, **settings).to_consistent_dict()
+ settings = self.DAQSettings(freq=freq, **settings)
+
+ if 'bandwidth' in settings:
+ warnings.warn('The bandwidth parameter has been replaced by f_max',
+ DeprecationWarning)
- if bandwidth == 'auto':
- bandwidth = settings['f_max'] / 2
+ if 'f_max' not in settings:
+ settings.f_max = settings.fs / 4
- if filter_order:
+ if filter_order is not None:
self.device.demods[self.demod].order(int(filter_order))
# BW 3dB = √(2^(1/n) - 1) / 2πτ
# BW NEP = Γ(n - 1/2) / 4τ √(π)Γ(n)
n = self.device.demods[self.demod].order()
- tc = gamma(n - 0.5) / (4 * bandwidth * np.sqrt(np.pi) * gamma(n))
+ tc = gamma(n - 0.5) / (4 * settings['f_max'] * np.sqrt(np.pi) * gamma(n))
# Do not use context manager here because somehow settings can get lost
# with device.set_transaction():