FFT vs least squares fitting of fourier components?

3
votes

So I've got a signal, and I've tried fitting a curve to it using two methods that I thought should have been numerically equivalent, but apparently are not.

Method 1: Explicit fitting of sinusoids by least squares:

def curve(x, a0, a1, b1, a2, b2):
  return a0 + a1*np.cos(x/720*2*math.pi) + b1*np.sin(x/720*2*math.pi) + a2*np.cos(x/720*2*math.pi*2) + b2*np.sin(x/720*2*math.pi*2)

def fit_curve(xdata, ydata):
  guess = [10, 0, 0, 0, 0]
  params, params_covariance = optimize.curve_fit(curve, xdata, ydata, guess)
  return params, params_covariance

Method 2: Use of inbuilt FFT algorithm to do the same thing:

  f = np.fft.rfft(y,3)
  curve = np.fft.irfft(f, width)

I have two problems. The first one is minor, the FFT is 'out of scale', so I apply a scaling factor mean(y)/mean(curve) to fix it, which is a bit of a hack. I'm not sure why this is the case.

The main problem I have is that I believe these should produce almost identical results, but they don't. The explicit fitting produces a tighter fit every time than the FFT results- my question is, should it?

2
pythonnumpyscipyfftcurve-fitting
I am currently consumed by this question, and what I have gathered till now tells me that least-squares fitting is somehow not utilising the real power of the Fourier transform, and hence less "strong". I can find no literature that establishes their equivalence. - Milind R

2 Answers

3
votes

Take a look at the docstring for np.fft.rfft. In particular, this: "If n is smaller than the length of the input, the input is cropped." When you do this:

    f = np.fft.rfft(y,3)

you are computing the FFT of the first three data points in y, not the first three Fourier coefficients of y.

3
votes

One can find discrete fourier transform coefficients using linear algebra, though I imagine it's only useful to understand the DFT better. The code below demonstrates this. To find coefficients and phases for a sine series will take a little bit more work, but shouldn't be too hard. The Wikipedia article cited in the code comments might help.

Note that one doesn't need scipy.optimize.curve_fit, or even linear least-squares. In fact, although I've used numpy.linalg.solve below, that is unnecessary since basis is a unitary matrix times a scale factor.

from __future__ import division, print_function

import numpy

# points in time series
n= 101
# final time (initial time is 0)
tfin= 10

# *end of changeable parameters*

# stepsize
dt= tfin/(n-1)
# sample count
s= numpy.arange(n)
# signal; somewhat arbitrary
y= numpy.sinc(dt*s)
# DFT
fy= numpy.fft.fft(y)
# frequency spectrum in rad/sample
wps= numpy.linspace(0,2*numpy.pi,n+1)[:-1]

# basis for DFT
# see, e.g., http://en.wikipedia.org/wiki/Discrete_Fourier_transform#equation_Eq.2
# and section "Properties -> Orthogonality"; the columns of 'basis' are the u_k vectors
# described there
basis= 1.0/n*numpy.exp(1.0j * wps * s[:,numpy.newaxis])

# reconstruct signal from DFT coeffs and basis
recon_y= numpy.dot(basis,fy)

# expect yerr to be "small"
yerr= numpy.max(numpy.abs(y-recon_y))
print('yerr:',yerr)

# find coefficients by fitting to basis
lin_fy= numpy.linalg.solve(basis,y)

# fyerr should also be "small"
fyerr= numpy.max(numpy.abs(fy-lin_fy))
print('fyerr',fyerr)

On my system this gives 

yerr: 2.20721480995e-14
fyerr 1.76885950227e-13

Tested on Ubuntu 14.04 with Python 2.7 and 3.4.