CODEWITHSUNDEEP
pythonscipytime-seriesfftmodel-fitting
Hooloovoo
Hooloovoo

Reputation: 875

Fitting time series with Fourier components: estimating Fourier series coefficients

Problem: I have a set of measurements (time, measurement, error) that exhibit periodic variations and I want to fit them with a Fourier series of the form

enter image description here

where A0 is the mean value of my measurements, t is the time, t0 is a (known) reference time and P is the (known) period. I want to fit for the coefficients A_k and phi_k.

Here is what I've got at the moment:

# Find Fourier components
# nfourier is the number of fourier components 

    def fourier(coeffs, time_data, epoch, period, nfourier, A0):
       import numpy as np
       omega = 2.0*np.pi/period
       fseries = np.zeros(len(time_data))
       fseries.fill(A0)
       for k in range(nfourier):
          ak = coeffs[k]
          phik = coeffs[k+1]
          time_diff = time_data - epoch
          fseries = fseries + ak * np.cos(k * omega * time_diff + phik)

       return fseries

I estimate the residuals as follows:

def residuals(coeffs, measurement_data, time_data, error_data, epoch, period, nfourier, A0):
   model = fourier(coeffs, time_data, epoch, period, nfourier, A0)
   result = measurement_data - model
   return result

Then I fit it with:

def fit_it(coeffs, measurement_data, time_data, error_data, epoch, period, nfourier, A0):
   from scipy.optimize import leastsq
   opt_coeff = leastsq(residuals, coeffs, args=(measurement_data, time_data, error_data, epoch, period, nfourier, A0))
   return opt_coeff

The program completes successfully but the fit seems to fail as can be seen from this figure: enter image description here

I am not sure what I am doing wrong here but maybe an expert can offer some advice. I would be happy to provide a test dataset if someone is willing to help.

Upvotes: 2

Views: 1892

Answers (2)

Michael Busch
Michael Busch

Reputation: 1

If you know the period of the data you should phase-fold your x-axis. It looks like the x-axis is in Julian Day. You should calculate the phase during your measurements.

Phase = ((Time - Reference Time) % Period) / Period

You will want to fit the fourier series to the measurement vs phase plot, which would look much more like a periodic signal.

Upvotes: 0

hsinghal
hsinghal

Reputation: 152

I dont understand the reason behind the Fourier fitting. I think you want to know the modulation frequency components of your data. I would suggest that you take the mean of data at each time and take fft of that, it will give you more insight about the frequency spectrum of your data.

As far as your code is concerned i would like to make two comments. First the phase of kth element is amplitude of k+1 th element. And second the error_data is not taking anything from function residuals. You may check over these points.

This is more like a comment but i dont have enough reputation to post comment. Just trying to help.

Regards

Upvotes: 1

Related Questions