CODEWITHSUNDEEP
pythonastropysymfit
Daryl Janzen
Daryl Janzen

Reputation: 1

Improved Fourier Series Parameters in Python

Edit: to clarify the problem, I am specifically looking for a routine that will accurately compute phi1, phi2, phi3, etc in the Fourier decomposition of periodic data, given a series of the form

m(t) = A0 + sum_{n=1}^{N} A_i * sin(2*pi*n*f + phi_n)

As discussed below, I've explored some of the related discussion on Stack Overflow, but none seem to produce accurate results (unless I've made an error in my code or math).

I'm trying to fit a Fourier series to some periodic variable star data, from which I want specific model parameters. Either equation (1) or (2) here, or anything equivalent would be useful (I apologise for using a link; I am not allowed to post pictures, and can't see how to typeset equations except as plain text). The results would be used, for instance, to calculate metallicity for an RRab Lyrae variable star, as described e.g. in this paper. Essentially, this requires an accurate determination of phi_31=phi_3-3*phi_1, in this case using the sine form of the Fourier series (equation 2).

There are numerous questions on this site about fitting Fourier series, but I don't think any of the answers quite get to what I need. The answer here actually uses symfit to model equation (1) from the Deb and Singh paper linked above. Using this method, I get a model that fits quite well to RRab Lyrae data, but the empirical results appear to be way off, indicating that the model doesn't find the right parameters. I've also modified that code to the sine Fourier series, and the best fit model yields different results, so I think symfit is just not doing a good enough job (my code is a bit long for this post, but I'd be happy to share it if anyone wants).

The answer given by jakevdp here indicates that symfit is getting stuck in a local minimum, which to me makes sense. The Lomb-Scargle approach he discusses does an excellent job, and since the code includes with his answer has some depricated functions I'll update it. First, we need some data to work with. Here's some photometric data for the RRab Lyrae variable SV Hya, from the Catalina Surveys DR2. The code below roughly follows jakevdp's code to pull in the data, plot the LS periodogram around the star's actual period (0.47855 days), and pull out the best-fit LS model parameters:

import numpy as np
data=np.genfromtxt('result_web_filePAevp8.csv', delimiter=',')
time=data[1:281,5]
signal=data[1:281,1]
signalerror=data[1:281,2]
import matplotlib.pyplot as plt
plt.figure()
from astropy.timeseries import LombScargle
ls = LombScargle(time, signal, signalerror, nterms=5)
freq, power = ls.autopower(minimum_frequency=2.05,maximum_frequency=2.15,samples_per_peak=10)
plt.plot(freq, power);
plt.show()
best_frequency = freq[np.argmax(power)]
theta=ls.model_parameters(best_frequency)
print(theta)

This is really cool, but my issue here is that I don't see how, as jakevdp says,

These linear sin/cos amplitudes can be converted back into the non-linear amplitude and phase with a bit of trigonometry.

I don't see how the Lomb-Scargle model, described here, can be rewritten as either equation (1) or (2) from Deb and Singh. The LS model has N+1 terms, not 2N+1, so I guess a 2*nterms LS model would be needed to yield an N-term Deb and Singh-type model. However, if you look at the model at the link, it's clear that setting cos(a)=sin(a+pi/2) enables us to rewrite the model in a form similar to equation (2) in Deb and Singh, only with the phi parameters fixed at either 0 for even terms and pi/2 for odd terms. Therefore, the model seems to constrain the actual parameters a priori that I want to determine.

So, my questions are these: (1) is my math wrong, and can I actually convert the LS model results to a form that gives me what I need? (2) If not, is there a way to edit the series in the LS model to the more general 2N+1-term form? (3) If not, is there another method I can use to fit the general 2N+1-term Fourier series to these data so I can extract phi1, phi3, etc?

Upvotes: 0

Views: 265

Answers (1)

Mike
Mike

Reputation: 1827

The question marked "Edit" is a well-known result from spectral analysis. Here is a straightforward routine that will compute phi1, phi2, etc. (Please be advised that the function below computes the "full" amplitude spectrum, so half the power is in negative frequencies, but the phase is what is relevant to your interests and is unaffected.)

import numpy as np
import matplotlib.pyplot as plt

def fft(t, x):
    fftspec = np.fft.fft(x)
    freqs = np.fft.fftfreq(x.size, np.mean(np.ediff1d(t)))
    amp = np.abs(fftspec)/len(fftspec)
    phase = np.arctan2(np.imag(fftspec), np.real(fftspec))
    #unfuck the phases to only show up
    #when amplitude is above some threshold
    phase = phase*(amp>1e-12)
    return freqs, amp, phase

if __name__ == "__main__":
    fs = 200.0 #Hz
    ampl1, ampl2   = 5, 7 #Volts
    freq1, freq2   = 15, 22 #Hz
    phase1, phase2 = 1, 2 #radians
    t = np.arange(0,10,1/fs) #timesteps, 10 seconds
    #2-component signal
    x1 = ampl1*np.cos(2*np.pi*freq1*t+phase1) #first component
    x2 = ampl2*np.cos(2*np.pi*freq2*t+phase2) #second component
    x = x1 + x2
    freqs, amp, phase = fft(t, x)
    fig, (ax0, ax1, ax2) = plt.subplots(nrows=1, ncols=3,
                                        figsize = (15,6))
    ax0.plot(t,x)
    ax0.set_title('Input Data')
    ax0.set_ylabel('Volts')
    ax0.set_xlabel('Time[s]')

    ax1.plot(freqs, amp)
    ax1.set_title('Full Amplitude Spectrum')
    ax1.set_ylabel('Volts')
    ax1.set_xlabel('Freq[Hz]')

    ax2.plot(freqs, phase)
    ax2.set_title('Phase spectrum')
    ax2.set_ylabel('Phase')
    ax2.set_xlabel('Freq[Hz]')

    plt.show()

Upvotes: 0

Related Questions