How to calculate amplitude from astropy Lomb-Scargle normalized psd?

Question

I'm using the Lomb-Scargle package in astropy. I tried with artificial data which is a sin function with amplitude of 1:

from astropy.stats import LombScargle
t_sin2=np.arange(1000)*1.0
a_sin2=np.sin(t_sin2)
frequency=np.arange(0.001,0.5,0.001)
PSD_LS = LombScargle(t_sin2,a_sin2).power(frequency, normalization='psd')
plt.plot(frequency, PSD_LS)

The plot I got is: my PSD plot. The peak value of PSD is around 230. I don't know how to calculate it to amplitude.

This is the usage of Lomb-Scargle in astropy:Lomb-Scargle docs. But I'm confused with the PSD normalization. In the usage, it says: explaination of PSD normalization, and χref is the best-fit sum-of-residuals of least-squares fits around a constant reference model, which is the term I do not understand.

Thank you!

Answer 1

power spectral density PSD is only the squared amplitude, hence if you calculate the square root of the PSD you can convert it to the amplitude:

import numpy as np    
Amp = np.sqrt(PSD_LS)

Answer 2

After digging into the original implementation of the code, it seems like PSD implementation takes into account amount of points you have as well.

    elif normalization == 'psd':
       p *= 0.5 * (dy ** -2.0).sum()

p are the calculated LS powers, while dy are the uncertainties for each data point. If you pass no uncertainties, it becomes an array of ones that has the same length as the amount of points that you have. In other words, it multiplies the power by the amount of points that you have (because of the .sum()). So to get an amplitude estimation, you take the PSD peak, divide by the amount of points, take square root and multiply by two. That should be an estimate of the peak.

sqrt(PSD_peak / amount_of_points) * 2 ? amplitude_of_the_fit

I am not sure why the astropy implementation multiplies by the sum of the uncertainties, because that results in this kind of behaviour when no uncertainties are given. Square root is needed because the units of the power are y.unit^2, just like in the original documentation.

For your example: sqrt(230 / 1000) * 2 = 0.959, which is close enough to the amplitude of 1. From my own trials on real data, I get a very close estimate of the amplitude fit for most objects, so that should hopefully be a correct formula.