cutcutcodec.core.signal.psd.welch

cutcutcodec.core.signal.psd.welch(signal_1: Tensor, signal_2: Tensor | None = None, freq_res: Real | None = None) → Tensor[source]

Estimate the power spectral density (PSD) ie intercorrelation with the Welch method.

Terminology

$x_1$ and $x_2$ are the signal whose we want to estimate the intercorrelation.
$s_x$ is the number of samples in the signal.
$s_w$ is the number of samples in the dpss window.
$r$ is the normalize frequency resolution in Hz, for a sample rate of 1.
$\Sigma(f)$ is an estimation of the standard deviation of the psd.
$\Gamma_{x_1,x_2}(f) = psd(f)$ is the power spectral density or the intercorrelation.
$n_{psd}(f)$ is the psd noise, ignoring gibbs effects.
$n_{win}$ is the additive noise created by the gibbs effect of the window.
$\eta$ is the maximum amplitude of the largest of the window’s secondary lobs.
$\alpha$ is the window theorical standardized half bandwidth.
$\beta$ is the window experimental standardized half bandwidth.

Equations

There, the equation to find the optimal windows size:

\[\begin{split}\begin{cases} n_{psd}(f) = \Sigma(f) . \sqrt{\frac{s_w}{s_x}} & \text{std of mean estimator} \\ n_{win} = \max\left(\Gamma_{x_1,x_2}(f)\right) . \eta & \text{because convolution}\\ \eta = g(\alpha) \\ \beta = \alpha + \epsilon = h(\alpha) \\ r = \frac{1}{s_w} + 2 . \frac{\beta}{s_w} & \text{convolution main lob} \\ \end{cases}\end{split}\]

To avoid having a too big gibbs noise, we want $n_{win} < n_{psd}(f)$.

\[\begin{split}\Leftrightarrow \begin{cases} psd_{max} . \eta < \Sigma_{max} . \sqrt{\frac{s_w}{s_x}} \\ r . s_w = 1 + 2.h(g^{-1}(\eta)) \\ \end{cases} \Leftrightarrow s_w . r - 2 . h\left( g^{-1}\left(\frac{\Sigma_{max} . \sqrt{\frac{s_w}{s_x}}}{psd_{max}}\right) \right) - 1 > 0\end{split}\]

Parameters

signal_1torch.Tensor: The stationary signal $x_1$ on witch we evaluate the PSD. The tensor can be batched, so the shape is (…, n).
signal_2torch.Tensor, default=signal_1: If you prefer to compute an intercorrelation rather an autocorrelation, please provide it ($x_2$), otherwise, $x_2 = x_1$.
freq_resfloat: The normlised frequency resolution in Hz $\left(r \in \left]0, \frac{1}{2}\right[\right)$, for a sample rate of 1. Higher it is, better is the frequency resolution but noiser it is.

Returns

psdtorch.Tensor: An estimation of the power spectral density $\Gamma_{x_1,x_2}(f)$, of shape (…, m). Ie an estimation of $\mathbb{E}\left[X_1(f)X_2^x(f)\right]$. In the case of autocorrelation, psd is returned as float type, overwise (in case of intercorrelation), it is returned as complex.

Notes

The complexity is $O\left(\frac{s_x}{r}\right)$.

Examples

>>> import torch
>>> from cutcutcodec.core.signal.psd import welch
>>> sr, delta_t = 8000, 60  # sample rate in (Hz) and duration in (s)
>>> t = torch.arange(0, delta_t, 1/sr)  # timae smaple in (s)
>>> signal = torch.randn((16, 2, len(t))) + torch.sin(2*torch.pi*440*t)  # 16 stereo signals
>>> psd = welch(signal, freq_res=20.0/sr)  # band of 20 Hz
>>>
>>> # we should observe a background at 1 and a pic a 440 Hz.
>>> # import matplotlib.pyplot as plt
>>> # freq = torch.fft.rfftfreq(2*psd.shape[-1]-1, 1/sr)
>>> # _ = plt.plot(freq, psd[0].T)
>>> # _ = plt.xlabel("freq (Hz)")
>>> # _ = plt.ylabel("psd")
>>> # plt.show()
>>>