python - 傅立叶分量的FFT与最小二乘拟合？

所以我得到了一个信号，我尝试使用两种我认为应该在数值上等效的方法来拟合曲线，但显然不是。

方法 1:通过最小二乘法显式拟合正弦曲线:

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

方法 2:使用内置 FFT 算法来做同样的事情 :

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

我有两个问题。第一个是次要的，FFT 是“超出规模”，所以我应用了一个缩放因子 mean(y)/mean(curve) 来修复它，这有点小技巧。我不确定为什么会这样。

我的主要问题是我相信这些应该产生几乎相同的结果，但他们没有。显式拟合每次都会产生比 FFT 结果更紧密的拟合——我的问题是， 应该是吗？

最佳答案

可以使用线性代数找到离散傅立叶变换系数，但我认为只有更好地理解 DFT 才有用。下面的代码演示了这一点。找到正弦序列的系数和相位需要做更多的工作，但应该不会太难。代码注释中引用的维基百科文章可能会有所帮助。

请注意，不需要 scipy.optimize.curve_fit ，甚至不需要线性最小二乘法。事实上，虽然我在下面使用了 numpy.linalg.solve，但这是不必要的，因为 basis 是一个酉矩阵乘以比例因子。

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)

在我的系统上，这给出了

yerr: 2.20721480995e-14
fyerr 1.76885950227e-13

在 Ubuntu 14.04 和 Python 2.7 和 3.4 上测试。

关于python - 傅立叶分量的FFT与最小二乘拟合？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/23561856/