2

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.

Hooloovoo
  • 865
  • 2
  • 11
  • 21
  • why dont use the definition and estimate the coefficients by averaging the scaled data with each `cos()` factor? should be simpler. Hard to understand the code you use – Nikos M. Dec 18 '15 at 17:55
  • Hi Niko, can you give me an example of what you mean? My data is just a time series and I'm passing some of the fixed parameters (epoch, period, nfourier, A0) in as arguments between the functions. The function is defined in fourier() and fit_it() minimizes (or at least it should) the residuals. It is just a linear least squares fit. – Hooloovoo Dec 18 '15 at 18:55
  • ok, i was refering the actual definition of [fourier coeficients](https://en.wikipedia.org/wiki/Fourier_series#Definition), but anyway, still do not understand the code too well, for example how is epoch related, or some other params? – Nikos M. Dec 18 '15 at 19:04
  • Based on the formula I gave in the original post: coeffs - These are the A_k, phi_k coefficients. measurement_data - These are my y-axis values. time_data - The x-axis values (t) error_data - The uncertainties for each point in measurement_data epoch - The known reference time (t0) period - The known period (P) nfourier - Number of fourier parameters to fit (e.g. 9) A0 - mean(measurement_data) Does this make sense? – Hooloovoo Dec 19 '15 at 00:47
  • 1
    using least squares as an approximation is an approach, however, note that least squares will try to find a global approximation. This means that for example by changing the `nfourier` param, for example, from 2 to 3 and doing least squares, the first 2 coefficients of `nfourier` 3 will NOT match the first 2 coefficients of `nfourier` 2 (while these should match). Again if i nunderstand the code correctly so far – Nikos M. Dec 19 '15 at 09:17

2 Answers2

1

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

hsinghal
  • 152
  • 1
  • 10
  • Hi hsinghal, Thanks for the comments. I will need to modify the code and see how it performs after I take them into account. As regards the first part of your comment, I need the fourier coefficients to calculate some other parameters. – Hooloovoo Dec 21 '15 at 15:57
0

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.

Michael Busch
  • 1
  • 1