Scipy periodogram terminology confusion

Question

I am confused about the terminology used in scipy.signal.periodogram, namely:

scaling : { 'density', 'spectrum' }, optional Selects between computing the power spectral density ('density') where Pxx has units of V*2/Hz if x is measured in V and computing the power spectrum ('spectrum') where Pxx has units of V*2 if x is measured in V. Defaults to 'density'

(see: http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.periodogram.html)

1) a few tests show that result for option 'density' is dependent on signal and window length and sampling frequency (grows when signal length increases). How come? I would say that it is exactly density that should be not dependent on these things. If I take a longer signal I should just get more accurate estimation, not different result. Not to mention that dependence on window length is also very surprising. Result diverges in the limit of infinite signal, which could be a feature of energy, but not power. Shouldn't the periodogram converge to real theoretical PSD when length increases? If, so, am I supposed to perform another normalisation outside of the signal.periodogram method?

2) to the contrary I see that alternative option 'spectrum' gives what I would previously call Power Spectrum Density, that is, it gives a resuls independent on window segment and window length and consistent with theoretical calculation. For instance for Asin(2PIft) a two sided solution yields two peaks at -f and f, each of height 0.25*A^2.

There is a lot of literature on this subject, but I get an impression that also there is a lot of incompatibile terminology, so I will be thankful for any clarification. The straightforward question is how to interpret these options and their units. (I am used to seeing V^2/Hz which are labeled "Power Spectrum Density").

Supranational

This video explains the difference youtu.be/… – Stagnate

Answer 1

Let's take a real array called data, of length N, and with sampling frequency fs. Let's call the time bin dt=1/fs, and T = N * dt. In frequency domain, the frequency bin df = 1/T = fs/N.

The power spectrum PS (scaling='spectrum' in scipy.periodogram) is calculated as follow:

import numpy as np
import scipy.fft as fft
dft = fft.fft(data)
PS = np.abs(dft)**2 / N ** 2

It has the units of V^2. It can be understood as follow. By analogy to the continuous Fourier transform, the energy E of the signal is:

E := np.sum(data**2) * dt = 1/N * np.sum(np.abs(dft)**2) * dt

(by Parseval's theorem). The power P of the signal is the total energy E divided by the duration of the signal T:

P := E/T = 1/N**2 * np.sum(np.abs(dft)**2)

The power P only depends on the Discrete Fourier Transform (DFT) and the number of samples N. Not directly on the sampling frequency fs or signal duration T. And the power per frequency channel, i.e., power spectrum SP, is thus given by the formula above:

PS = np.abs(dft)**2 / N ** 2

For the power spectrum density PSD (scaling='density' in scipy.periodogram), one needs to divide the PS by the frequency bin of the DFT, df:

PSD := PS/df = PS * N * dt = PS * N / fs

and thus:

PSD = np.abs(dft)**2 / N * dt

This has the units of V^2/Hz = V^2 * s, and now depends on the sampling frequency. That way, integrating the PSD over the frequency range gives the same result as summing the individual values of the PS.

This should explain the relations that you see when changing the window, sampling frequency, duration.

Padishah

Answer 2

1. scipy.signal.peridogram uses the scipy.signal.welch function with 0 overlap. Therefore, the scaling is similar to the one provided by the welch function, density or spectrum.
2. In case of the density scaling, the amplitude will vary with window length, as the longer the window the higher the frequency resolution e.g. the \Delta_f is smaller. Since the estimated density is the average one, the smaller the \Delta_f the less zero energy is considered in the averaging.
3. As you have mentioned spectrum scaling is an integration of the energy density over the spectrum to produce the energy. Therefore, the integration over zero values does not affect the final value.

Scenic

Answer 3

1. Fourier transform actually requires finite energy in an infinite duration of time series (like a decay). So, If you just make your time series sample longer by "duplicating", the energy will be infinite with an infinite duration.
2. My main confusion was on the "spectrum" option for scipy.signal.periodogram, which seems to create a constant energy spectrum even when the time series become longer. Normally, 0.5*A^2=S(f)*delta_f, where S(f) is the power density spectrum. S(f)*delta_f, representing energy is constant if A is constant. But when using a longer duration of time series, delta_f (i.e. incremental frequency) is reduced accordingly, based on FFT procedure. For example, 100s time series will lead to a delta_f=0.01Hz, while 1000s time series will have a delta_f=0.001Hz. S(f) representing density will accordingly change.

Variola