cutcutcodec.core.signal.psd.welch

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

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

Terminology

  • \(s_x\) is the number of samples in the signal.

  • \(s_w\) is the number of samples in the dpss window.

  • \(r\) is the frequency resolution in Hz, for a sample rate of 1.

  • \(\sigma(f)\) is an estimation of the std of the psd.

  • \(psd(f)\) is the power spectral density.

  • \(n_{psd}(f)\) is the psd noise, ignoring gibbs effects.

  • \(n_{win}\) is the additive noise created by the gibbs effect of the window.

  • \(\gamma\) 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(psd(f)\right) . \gamma & \text{because convolution}\\ \gamma = 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} . \gamma < \sigma_{max} . \sqrt{\frac{s_w}{s_x}} \\ r . s_w = 1 + 2.h(g^{-1}(\gamma)) \\ \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

signaltorch.Tensor

The stationary signal on witch we evaluate the PSD. The tensor can be batched, so the shape is (…, n).

freq_resfloat

The normlised frequency resolution in Hz (\(r \in \left]0, \frac{1}{2}\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, of shape (…, m).

Examples

>>> import torch
>>> from cutcutcodec.core.signal.psd import welch
>>> signal = torch.randn((16, 2, 768000))  # 16 stereos signals
>>> psd = welch(signal, 1/4800)
>>>
>>> # freq = torch.fft.rfftfreq(2*psd.shape[-1]-1, 1/48000)
>>> # import matplotlib.pyplot as plt
>>> # _ = plt.plot(freq, psd[0].T)
>>> # plt.show()
>>>