假设我有一些数据 y，我想对其进行傅立叶级数拟合。对此post https://stackoverflow.com/questions/4258106/how-to-calculate-a-fourier-series-in-numpy，解决方案由Mermoz使用级数的复杂格式并“用黎曼和计算系数”。在这另一post https://dsp.stackexchange.com/questions/40780/how-to-get-the-fourier-series-using-pythons-tt-fft/49233，通过FFT得到级数，并写出一个例子。
我尝试实现这两种方法（下面的图像和代码 - 请注意，每次运行代码时，由于使用numpy.随机.正常）但我想知道为什么我得到不同的结果 - 黎曼方法似乎“错误地转移”，而 FFT 方法似乎“挤压”。我也不确定我对该系列的周期“tau”的定义。我很感谢您的关注。
我在 Windows 7 上使用 Spyder 和 Python 3.7.1
Example https://i.stack.imgur.com/bjJ6T.png
import matplotlib.pyplot as plt
import numpy as np
# Assume x (independent variable) and y are the data.
# Arbitrary numerical values for question purposes:
start = 0
stop = 4
mean = 1
sigma = 2
N = 200
terms = 30 # number of terms for the Fourier series
x = np.linspace(start,stop,N,endpoint=True)
y = np.random.normal(mean, sigma, len(x))
# Fourier series
tau = (max(x)-min(x)) # assume that signal length = 1 period (tau)
# From ref 1
def cn(n):
c = y*np.exp(-1j*2*n*np.pi*x/tau)
return c.sum()/c.size
def f(x, Nh):
f = np.array([2*cn(i)*np.exp(1j*2*i*np.pi*x/tau) for i in range(1,Nh+1)])
return f.sum()
y_Fourier_1 = np.array([f(t,terms).real for t in x])
# From ref 2
Y = np.fft.fft(y)
np.put(Y, range(terms+1, len(y)), 0.0) # zero-ing coefficients above "terms"
y_Fourier_2 = np.fft.ifft(Y)
# Visualization
f, ax = plt.subplots()
ax.plot(x,y, color='lightblue', label = 'artificial data')
ax.plot(x, y_Fourier_1, label = ("'Riemann' series fit (%d terms)" % terms))
ax.plot(x,y_Fourier_2, label = ("'FFT' series fit (%d terms)" % terms))
ax.grid(True, color='dimgray', linestyle='--', linewidth=0.5)
ax.set_axisbelow(True)
ax.set_ylabel('y')
ax.set_xlabel('x')
ax.legend()
执行两个小的修改足以使总和与 np.fft 的输出几乎相似。FFTW 库确实计算了这些总和 http://www.fftw.org/fftw3_doc/The-1d-Discrete-Fourier-Transform-_0028DFT_0029.html#The-1d-Discrete-Fourier-Transform-_0028DFT_0029.
1) 信号的平均值，c[0]应计入：
f = np.array([2*cn(i)*np.exp(1j*2*i*np.pi*x/tau) for i in range(0,Nh+1)]) # here : 0, not 1
2) 输出必须按比例缩放。
y_Fourier_1=y_Fourier_1*0.5
The output seems "squeezed" because the high frequency components have been filtered. Indeed, the high frequency oscillations of the input have been cleared and the output looks like a moving average.
Here, tau实际上定义为stop-start：它对应于帧的长度。这是信号的预期周期。
如果帧与信号的周期不对应，您可以通过将信号与其自身进行卷积并找到第一个最大值来猜测其周期。看从 FFT 中查找信号的周期 https://stackoverflow.com/questions/49531952/find-period-of-a-signal-out-of-the-fft/49545394#49545394然而，它不太可能与由生成的数据集正常工作numpy.random.normal： 这是一加性高斯白噪声 https://en.wikipedia.org/wiki/Additive_white_Gaussian_noise。由于它具有恒定的功率谱密度，因此很难将其描述为周期性的！