Bug with spectra when using real fft

[dougiesquire]

I think there is a bug with the power spectrum estimates when a real dimension is provided.

For real fft, Hermitian symmetry is exploited to compute only the non-negative frequency terms. But we can't ignore the energy in the negative frequency terms when computing spectra. To account for this, we need to multiply by 2 the energy of all but the first (DC) mode (and last mode if the signal length even - see https://numpy.org/doc/stable/reference/generated/numpy.fft.rfft.html).

I can demonstrate the issue and the fix through comparison to scipy.signal.periodogram:

import xrft
import numpy as np
import scipy.signal
import xarray as xr
import matplotlib.pyplot as plt 

x = xr.DataArray(np.random.normal(0, 1, size=100), 
                 dims=['t'], 
                 coords=[range(0,100)])


f_scipy, p_scipy = scipy.signal.periodogram(x.values,
                                            window='rectangular')
p_xrft = xrft.xrft.power_spectrum(x, 
                                  dim=['t'], 
                                  real='t',
                                  detrend='constant')

# For even length signals the last mode includes both the +ve and -ve contributions
correction = np.ones_like(f_scipy)
correction[1:-1] = 2

plt.plot(f_scipy, p_scipy,
         label='scipy psd')
plt.plot(p_xrft.freq_t, p_xrft, 
         color='C2', label='xrft psd')
plt.plot(p_xrft.freq_t, correction*p_xrft, 
         linestyle='--', color='C1', label='corrected xrft psd')
plt.xlabel('freq')
plt.ylabel('psd')
plt.legend()

[roxyboy]

These are all really great fixes! Thanks for putting in the time.

[lanougue]

Hi @dougiesquire, @roxyboy

I think this is not an xrft nor scipy periodogram bug. In scipy.signal.periodogram(), there is a default parameter "return_onesided" which is set to "True" as default. xrft power_spectrum() returns a two_sided spectrum.
If you set return_onesided parameter to False in scipy.signal.periodogram(), you should get a two_sided spectrum whose amplitude match xrft.power_spectrum.
Having return_onesided=True in scipy periodrogram() double the amplitude to compensate the "one sided".

[lanougue]

Actually, adding a "real" parameter in xrft makes xrft.power_spectrum() to return a one_sided spectrum whose amplitude is NOT doubled. real parameter in xrft is thus NOT an equivalent to return_onesided=True in scipy.periodogram() even if the two methods will return a one_sided spectrum.
There is thus no bug, it is just a question of convention that we maybe have to discuss.

[rabernat]

Thanks for the explanations @lanougue! I think it's super important we get all of this documented in the documentation.

[dougiesquire]

Thanks for your input, @lanougue.

I agree that it was probably not fair for me to call this a "bug". But, I think most users would expect that the output of a function for computing power spectrum should integrate to (approximately) the signal energy. This is the behaviour of scipy.signal.periodogram, for example, regardless of whether return_onesided is True or False. However, running xrft.power_spectrum with a real dimension currently returns a one-sided power spectrum that has not had its energy compensated and so integrates to only half of the signal energy. I personally think this is undesirable default behaviour.

Regardless of what is decided here, I agree with @rabernat that clear documentation is what is really important.

[lanougue]

Hi @dougiesquire,
I fully agree with you. We have to decide if we want to be compliant with scipy.signal.periodogram(). And populate the documentation for sure.