Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Inconsistent treatment of DC component #33

Open
keatonb opened this issue Oct 9, 2018 · 2 comments
Open

Inconsistent treatment of DC component #33

keatonb opened this issue Oct 9, 2018 · 2 comments

Comments

@keatonb
Copy link

keatonb commented Oct 9, 2018

I am attempting to compute a spectral window by taking the power spectrum of ones sampled as the time series. I expect a value of 1 at frequency=0, but usually I get zero power at zero frequency... but not always! I set fit_offset=False and center_data=False to try to retain the the DC value. Here's an example:

import numpy as np
import matplotlib.pyplot as plt
import gatspy.periodic as gp

time = np.linspace(0,9,1e5)
nyq=1./(2.*np.median(np.diff(time)))
deltaf=(1./time[-1])/4.5
freq,sw=gp.lomb_scargle_fast.lomb_scargle_fast(time,np.ones(time.size)/time.size,
                                               f0=0.,df=deltaf,Nf=nyq/deltaf,fit_offset=False,
                                               center_data=False)

plt.plot(freq,sw,label='oversample by 4.5')

deltaf=(1./time[-1])/5.
freq,sw=gp.lomb_scargle_fast.lomb_scargle_fast(time,np.ones(time.size)/time.size,
                                               f0=0.,df=deltaf,Nf=nyq/deltaf,fit_offset=False,
                                               center_data=False)
plt.plot(freq,sw,label='oversample by 5')
plt.title('why this different behavior at zero frequency?')
plt.legend()
plt.xlim(0,.5)
plt.xlabel('frequency')
plt.ylabel('spectral window')
plt.show()

gatspy_weirddc

This seems to only affect the power at exactly frequency = 0, and the result appears to be sensitive to some precise relationship between the frequency sampling and the time sampling. Am I correct to expect power=1 for frequency=0? Is the "spectral window" a well defined quantity for the Lomb-Scargle periodogram (though I expect it to reproduce the Fourier transform result for the evenly-sampled example above)?

@keatonb keatonb closed this as completed Oct 9, 2018
@keatonb keatonb reopened this Oct 9, 2018
@keatonb
Copy link
Author

keatonb commented Oct 9, 2018

Perhaps this is a result of the hardcoded

if np.isnan(power[0]) or np.isinf(power[0]):
    power[0] = 0

where that condition only happens sometimes. Should the forced behavior be different when center_data=False?

@keatonb
Copy link
Author

keatonb commented Oct 9, 2018

I believe the desired power at f = 0 is
power[0] = np.mean(y ** 2 - (y - np.mean(y)) ** 2) / np.mean(y ** 2)

I'll open a pull request.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant