从时间序列数据中确定傅立叶系数

Question

我问了一个关于如何从时间序列数据中确定傅里叶系数的已删除问题。我重新提交此文件是因为我已经更好地阐述了问题并有一个我会给出的解决方案，因为我认为其他人可能会发现这非常有用。

我有一些时间序列数据，我已将它们分箱到等间距的时间箱中（这对我的解决方案至关重要），并且我想从这些数据中确定最好的傅立叶级数（或任何函数）描述数据。这是一个带有一些测试数据的 MWE，以显示我要拟合的数据：

import numpy as np
import matplotlib.pyplot as plt

# Create a dependent test variable to define the x-axis of the test data.
test_array = np.linspace(0, 1, 101) - 0.5

# Define some test data to try to apply a Fourier series to.
test_data = [0.9783883464566918, 0.979599093567252, 0.9821424606299206, 0.9857575507812502, 0.9899278899999995,
             0.9941848228346452, 0.9978438300395263, 1.0003009205426352, 1.0012208923679058, 1.0017130521235522,
             1.0021799664031628, 1.0027475606936413, 1.0034168260869563, 1.0040914266144825, 1.0047781181102355,
             1.005520348837209, 1.0061899214145387, 1.006846206627681, 1.0074483048543692, 1.0078691461988312,
             1.008318736328125, 1.008446947572815, 1.00862051262136, 1.0085134881422921, 1.008337095516569,
             1.0079539881889774, 1.0074857334630352, 1.006747783037474, 1.005962048923679, 1.0049115434782612,
             1.003812267822736, 1.0026427549407106, 1.001251963531669, 0.999898555335968, 0.9984976286266923,
             0.996995982142858, 0.9955652088974847, 0.9941647321428578, 0.9927727076023389, 0.9914750532544377,
             0.990212467710371, 0.9891098035363466, 0.9875998927875242, 0.9828093773946361, 0.9722532524271845,
             0.9574084365384614, 0.9411012303149601, 0.9251820309477757, 0.9121488392156851, 0.9033119748549322,
             0.9002445803921568, 0.9032760564202343, 0.91192435882353, 0.9249696964980555, 0.94071381372549,
             0.957139088974855, 0.9721083392156871, 0.982955287937743, 0.9880613320235758, 0.9897455322896282,
             0.9909590626223097, 0.9922601592233015, 0.9936513112840472, 0.9951442427184468, 0.9967071285988475,
             0.9982921493123781, 0.9998775465116277, 1.001389230174081, 1.0029109110251453, 1.0044033691406251,
             1.0057110841487276, 1.0069551867704276, 1.008118776264591, 1.0089884470588228, 1.0098663972602735,
             1.0104514566473979, 1.0109849223300964, 1.0112043902912626, 1.0114717968750002, 1.0113343036750482,
             1.0112205972495087, 1.0108811786407768, 1.010500276264591, 1.0099054552529192, 1.009353759223301,
             1.008592596116505, 1.007887223091976, 1.0070715634615386, 1.0063525891472884, 1.0055587861271678,
             1.0048733732809436, 1.0041832862669238, 1.0035913326848247, 1.0025318871595328, 1.000088536345776,
             0.9963596140350871, 0.9918380684931506, 0.9873937281553398, 0.9833394624277463, 0.9803621496062999,
             0.9786476100386117]

# Create a figure to view the data.
fig, ax = plt.subplots(1, 1, figsize=(6, 6))

# Plot the data.
ax.scatter(test_array, test_data, color="k", s=1)

这会输出以下内容：

问题是如何确定最能描述此数据的傅里叶级数。确定傅里叶系数的常用公式需要在积分中插入一个函数，但如果我有一个函数来描述数据，我根本不需要傅里叶系数；找到这个系列的重点是要有数据的功能表示。那么，在没有这样一个函数的情况下，系数是怎么求出来的呢？

Answer 1

我对这个问题的解决方案是使用 NumPy 的快速傅里叶变换实现 numpy.fft.fft() 对数据应用离散傅里叶变换；这就是为什么数据在时间上均匀分布很重要的原因，因为 FFT 需要这样做。虽然 FFT 通常用于执行频谱分析，但所需的傅里叶系数与此函数的输出直接相关。

具体来说，这个函数输出一系列i复值系数c。傅立叶级数系数​​是使用以下关系求得的：

因此，FFT 允许直接计算傅立叶系数。这是我对这个问题的解决方案的 MWE，扩展了上面给出的示例：

import numpy as np
import matplotlib.pyplot as plt

# Set the number of equal-time bins to create.
n_bins = 101

# Set the number of Fourier coefficients to use.
n_coeff = 51

# Define a function to generate a Fourier series based on the coefficients determined by the Fast Fourier Transform.
# This also includes a series of phases x to pass through the function.
def create_fourier_series(x, coefficients):

    # Begin the series with the zeroeth-order Fourier coefficient.
    fourier_series = coefficients[0][0] / 2

    # Now generate the first through n_coeff'th terms.  The period is defined to be 1 since we're operating in phase
    # space.
    for n in range(1, n_coeff):
        fourier_series += (fourier_coeff[n][0] * np.cos(2 * np.pi * n * x) + fourier_coeff[n][1] * 
                           np.sin(2 * np.pi * n * x))

    return fourier_series

# Create a dependent test variable to define the x-axis of the test data.
test_array = np.linspace(0, 1, n_bins) - 0.5

# Define some test data to try to apply a Fourier series to.
test_data = [0.9783883464566918, 0.979599093567252, 0.9821424606299206, 0.9857575507812502, 0.9899278899999995,
             0.9941848228346452, 0.9978438300395263, 1.0003009205426352, 1.0012208923679058, 1.0017130521235522,
             1.0021799664031628, 1.0027475606936413, 1.0034168260869563, 1.0040914266144825, 1.0047781181102355,
             1.005520348837209, 1.0061899214145387, 1.006846206627681, 1.0074483048543692, 1.0078691461988312,
             1.008318736328125, 1.008446947572815, 1.00862051262136, 1.0085134881422921, 1.008337095516569,
             1.0079539881889774, 1.0074857334630352, 1.006747783037474, 1.005962048923679, 1.0049115434782612,
             1.003812267822736, 1.0026427549407106, 1.001251963531669, 0.999898555335968, 0.9984976286266923,
             0.996995982142858, 0.9955652088974847, 0.9941647321428578, 0.9927727076023389, 0.9914750532544377,
             0.990212467710371, 0.9891098035363466, 0.9875998927875242, 0.9828093773946361, 0.9722532524271845,
             0.9574084365384614, 0.9411012303149601, 0.9251820309477757, 0.9121488392156851, 0.9033119748549322,
             0.9002445803921568, 0.9032760564202343, 0.91192435882353, 0.9249696964980555, 0.94071381372549,
             0.957139088974855, 0.9721083392156871, 0.982955287937743, 0.9880613320235758, 0.9897455322896282,
             0.9909590626223097, 0.9922601592233015, 0.9936513112840472, 0.9951442427184468, 0.9967071285988475,
             0.9982921493123781, 0.9998775465116277, 1.001389230174081, 1.0029109110251453, 1.0044033691406251,
             1.0057110841487276, 1.0069551867704276, 1.008118776264591, 1.0089884470588228, 1.0098663972602735,
             1.0104514566473979, 1.0109849223300964, 1.0112043902912626, 1.0114717968750002, 1.0113343036750482,
             1.0112205972495087, 1.0108811786407768, 1.010500276264591, 1.0099054552529192, 1.009353759223301,
             1.008592596116505, 1.007887223091976, 1.0070715634615386, 1.0063525891472884, 1.0055587861271678,
             1.0048733732809436, 1.0041832862669238, 1.0035913326848247, 1.0025318871595328, 1.000088536345776,
             0.9963596140350871, 0.9918380684931506, 0.9873937281553398, 0.9833394624277463, 0.9803621496062999,
             0.9786476100386117]

# Determine the fast Fourier transform for this test data.
fast_fourier_transform = np.fft.fft(test_data[n_bins / 2:] + test_data[:n_bins / 2])

# Create an empty list to hold the values of the Fourier coefficients.
fourier_coeff = []

# Loop through the FFT and pick out the a and b coefficients, which are the real and imaginary parts of the
# coefficients calculated by the FFT.
for n in range(0, n_coeff):
    a = 2 * fast_fourier_transform[n].real / n_bins
    b = -2 * fast_fourier_transform[n].imag / n_bins
    fourier_coeff.append([a, b])

# Create the Fourier series approximating this data.
fourier_series = create_fourier_series(test_array, fourier_coeff)

# Create a figure to view the data.
fig, ax = plt.subplots(1, 1, figsize=(6, 6))

# Plot the data.
ax.scatter(test_array, test_data, color="k", s=1)

# Plot the Fourier series approximation.
ax.plot(test_array, fourier_series, color="b", lw=0.5)

这会输出以下内容：

请注意，我如何定义 FFT（导入数据的后半部分，然后导入前半部分）是该数据生成方式的结果。具体来说，数据从 -0.5 运行到 0.5，但 FFT 假设它从 0.0 运行到 1.0，因此需要进行这种转变。

我发现这对于不包含非常尖锐和狭窄不连续性的数据非常有效。我很想知道是否有人对此问题有其他建议的解决方案，我希望人们发现这个解释清楚且有帮助。

Answer 2

不确定它是否对您有帮助；我写了一个程序来插入你的数据。这是使用 buildingblocks==0.0.15

完成的

请看下文，

import matplotlib.pyplot as plt
from buildingblocks import bb
import numpy as np

Ydata = [0.9783883464566918, 0.979599093567252, 0.9821424606299206, 0.9857575507812502, 0.9899278899999995,
             0.9941848228346452, 0.9978438300395263, 1.0003009205426352, 1.0012208923679058, 1.0017130521235522,
             1.0021799664031628, 1.0027475606936413, 1.0034168260869563, 1.0040914266144825, 1.0047781181102355,
             1.005520348837209, 1.0061899214145387, 1.006846206627681, 1.0074483048543692, 1.0078691461988312,
             1.008318736328125, 1.008446947572815, 1.00862051262136, 1.0085134881422921, 1.008337095516569,
             1.0079539881889774, 1.0074857334630352, 1.006747783037474, 1.005962048923679, 1.0049115434782612,
             1.003812267822736, 1.0026427549407106, 1.001251963531669, 0.999898555335968, 0.9984976286266923,
             0.996995982142858, 0.9955652088974847, 0.9941647321428578, 0.9927727076023389, 0.9914750532544377,
             0.990212467710371, 0.9891098035363466, 0.9875998927875242, 0.9828093773946361, 0.9722532524271845,
             0.9574084365384614, 0.9411012303149601, 0.9251820309477757, 0.9121488392156851, 0.9033119748549322,
             0.9002445803921568, 0.9032760564202343, 0.91192435882353, 0.9249696964980555, 0.94071381372549,
             0.957139088974855, 0.9721083392156871, 0.982955287937743, 0.9880613320235758, 0.9897455322896282,
             0.9909590626223097, 0.9922601592233015, 0.9936513112840472, 0.9951442427184468, 0.9967071285988475,
             0.9982921493123781, 0.9998775465116277, 1.001389230174081, 1.0029109110251453, 1.0044033691406251,
             1.0057110841487276, 1.0069551867704276, 1.008118776264591, 1.0089884470588228, 1.0098663972602735,
             1.0104514566473979, 1.0109849223300964, 1.0112043902912626, 1.0114717968750002, 1.0113343036750482,
             1.0112205972495087, 1.0108811786407768, 1.010500276264591, 1.0099054552529192, 1.009353759223301,
             1.008592596116505, 1.007887223091976, 1.0070715634615386, 1.0063525891472884, 1.0055587861271678,
             1.0048733732809436, 1.0041832862669238, 1.0035913326848247, 1.0025318871595328, 1.000088536345776,
             0.9963596140350871, 0.9918380684931506, 0.9873937281553398, 0.9833394624277463, 0.9803621496062999,
             0.9786476100386117]


Xdata=list(range(0,len(Ydata)))
Xnew=list(np.linspace(0,len(Ydata),200))
Ynew=bb.interpolate(Xdata,Ydata,Xnew,40)

plt.figure()
plt.plot(Xdata,Ydata)
plt.plot(Xnew,Ynew,'*')
plt.legend(['Given Data', 'Interpolated Data'])

plt.show()

如果你想进一步写代码，我也给出了代码，你可以看到源代码并学习：

import module
import inspect
src = inspect.getsource(module)
print(src)