FFT与傅立叶分量的最小二乘拟合？答案

【问题标题】：FFT vs least squares fitting of fourier components?FFT与傅立叶分量的最小二乘拟合？
【发布时间】：2014-06-27 00:44:33
【问题描述】：

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

方法一：用最小二乘法显式拟合正弦曲线：

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 结果更紧密的拟合 - 我的问题是，应该吗？

【问题讨论】：

我目前正被这个问题所困扰，到目前为止我收集到的信息告诉我，最小二乘拟合不知何故没有利用傅里叶变换的真正力量，因此不太“强”。我找不到任何文献证明它们是等价的。

标签： python numpy scipy fft curve-fitting

【解决方案1】：

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

请注意，不需要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 进行测试。

【讨论】：

【解决方案2】：

看看docstring for np.fft.rfft。特别是：“如果 n 小于输入的长度，则裁剪输入。”当你这样做时：

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

您正在计算 y 中前三个数据点的 FFT，而不是 y 的前三个傅立叶系数。

【讨论】：