# -*- coding: utf-8 -*-
"""
This file contains qudi methods for mathematical operations/transformations.
.. Copyright (c) 2021, the qudi developers. See the AUTHORS.md file at the top-level directory of this
.. distribution and on <https://github.com/Ulm-IQO/qudi-core/>
..
.. This file is part of qudi.
..
.. Qudi is free software: you can redistribute it and/or modify it under the terms of
.. the GNU Lesser General Public License as published by the Free Software Foundation,
.. either version 3 of the License, or (at your option) any later version.
..
.. Qudi is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
.. without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
.. See the GNU Lesser General Public License for more details.
..
.. You should have received a copy of the GNU Lesser General Public License along with qudi.
.. If not, see <https://www.gnu.org/licenses/>.
"""
__all__ = ("compute_ft", "ft_windows")
import numpy as np
from scipy.signal import windows as window_func
# Available windows to be applied on signal data before FT.
# To find out the amplitude normalization factor check either the scipy implementation on
# https://github.com/scipy/scipy/blob/v0.15.1/scipy/signal/windows.py#L336
# or just perform a sum of the window (oscillating parts of the window should be averaged out and
# constant offset factor will remain):
# MM=1000000 # choose a big number
# print(sum(signal.hanning(MM))/MM)
ft_windows = {'none': {'func': np.ones, 'ampl_norm': 1.0},
'hamming': {'func': window_func.hamming, 'ampl_norm': 1.0/0.54},
'hann': {'func': window_func.hann, 'ampl_norm': 1.0/0.5},
'blackman': {'func': window_func.blackman, 'ampl_norm': 1.0/0.42},
'triang': {'func': window_func.triang, 'ampl_norm': 1.0/0.5},
'flattop': {'func': window_func.flattop, 'ampl_norm': 1.0/0.2156},
'bartlett': {'func': window_func.bartlett, 'ampl_norm': 1.0/0.5},
'parzen': {'func': window_func.parzen, 'ampl_norm': 1.0/0.375},
'bohman': {'func': window_func.bohman, 'ampl_norm': 1.0/0.4052847},
'blackmanharris': {'func': window_func.blackmanharris, 'ampl_norm': 1.0/0.35875},
'nuttall': {'func': window_func.nuttall, 'ampl_norm': 1.0/0.3635819},
'barthann': {'func': window_func.barthann, 'ampl_norm': 1.0/0.5}
}
[docs]
def compute_ft(x_val, y_val, zeropad_num=0, window='none', base_corr=True, psd=False):
"""
Compute the Discrete Fourier Transform (DFT) or the Power Spectral Density (PSD)
of the power spectral density.
Parameters
----------
x_val : numpy.array
1D array representing the x-values.
y_val : numpy.array
1D array of the same size as `x_val`, representing the y-values.
zeropad_num : int, optional
Zero-padding (adding zeros to the end of the array). If `zeropad_num >= 0`,
the size of the array is extended by adding zeros to the end of `y_val`
before performing the Fourier Transform (FT). The resulting array will have
the length `(len(y_val)/2)*(zeropad_num+1)`. Zero-padding does not add more
information to the DFT; it only interpolates between the DFT values. Set
`zeropad_num=1` to obtain output arrays of the same size as the input arrays.
Default is `zeropad_num=0`.
window : str, optional
The window function to be applied to the y-values before calculating the FT.
If not specified, no window function is applied.
base_corr : bool, optional
Whether baseline correction should be performed before calculating the FT.
Default is `False`.
psd : bool, optional
Whether to compute the Power Spectral Density (PSD) instead of the Discrete
Fourier Transform (DFT). The PSD is the FT of the absolute square of the
y-values. Default is `False`.
Returns
-------
tuple of numpy.array
A tuple containing the DFT or PSD results:
- dft_x : numpy.array
The x-values of the DFT/PSD.
- dft_y : numpy.array
The y-values of the DFT/PSD.
Note that the length of the return arrays depends on `zeropad_num`:
`len(dft_x) = len(dft_y) = (len(y_val)/2)*(zeropad_num+1)`.
Notes
-----
- The return values of the FT have only half of the entries compared to the used signal input
if `zeropad_num=0`.
- In general, a window function should be applied to the y-data before calculating the FT to
reduce spectral leakage. The Hann window is often a good choice, e.g., `window='hann'`.
- Keep in mind the relation to the Fourier transform space:
T = delta_t * N_samples
where `delta_t` is the distance between time points and `N_samples` is the number of points
in the time domain. Consequently, the sample rate is:
f_samplerate = T / N_samples
The FT returns values from 0 to `f_samplerate`, or equivalently `-f_samplerate/2` to
`f_samplerate/2`.
- Difference between PSD and DFT:
The Power Spectral Density (PSD) describes how the power of your signal is distributed
over frequency, while the DFT shows the spectral content of your signal, i.e., the
amplitude and phase of harmonics in your signal.
"""
x_val = np.array(x_val)
y_val = np.array(y_val)
# Make a baseline correction to avoid a constant offset near zero frequencies. Offset of the
# y_val from mean corresponds to half the value at fft_y[0]
corrected_y = y_val
if base_corr:
corrected_y = y_val - y_val.mean()
ampl_norm_fact = 1.0
# apply window to data to account for spectral leakage:
if window in ft_windows:
window_val = ft_windows[window]["func"](len(y_val))
corrected_y = corrected_y * window_val
# to get the correct amplitude in the amplitude spectrum
ampl_norm_fact = ft_windows[window]["ampl_norm"]
# zeropad for sinc interpolation:
zeropad_arr = np.zeros(len(corrected_y) * (zeropad_num + 1))
zeropad_arr[: len(corrected_y)] = corrected_y
# Get the amplitude values from the fourier transformed y values.
fft_y = np.abs(np.fft.fft(zeropad_arr))
# Power spectral density (PSD) or just amplitude spectrum of fourier signal:
power_value = 1.0
if psd:
power_value = 2.0
# The factor 2 accounts for the fact that just the half of the spectrum was taken. The
# ampl_norm_fact is the normalization factor due to the applied window function (the offset
# value in the window function):
fft_y = ((2 / len(y_val)) * fft_y * ampl_norm_fact) ** power_value
# Due to the sampling theorem you can only identify frequencies at half of the sample rate,
# therefore the FT contains an almost symmetric spectrum (the asymmetry results from aliasing
# effects). Therefore take the half of the values for the display.
middle = int((len(zeropad_arr) + 1) // 2)
# sample spacing of x_axis, if x is a time axis than it corresponds to a timestep:
x_spacing = np.round(x_val[-1] - x_val[-2], 12)
# use the helper function of numpy to calculate the x_values for the fourier space. That
# function will handle an occuring devision by 0:
fft_x = np.fft.fftfreq(len(zeropad_arr), d=x_spacing)
return abs(fft_x[:middle]), fft_y[:middle]
def normalize(arr: np.ndarray, axis=-1, order=2) -> np.ndarray:
"""
Taken from stack overflow
https://stackoverflow.com/questions/21030391/how-to-normalize-a-numpy-array-to-a-unit-vector
"""
l2 = np.atleast_1d(np.linalg.norm(arr, order, axis))
l2[l2 == 0] = 1
return arr / np.expand_dims(l2, axis)