python - 计算傅里叶级数来表示数据点

Question

我希望计算一个通过一组给定点的函数（傅立叶级数）。

类似于这里发生的事情https://gofigure.impara.ai/，但我不希望对其进行动画处理。我只想要这个功能，以便我可以自己绘制形状。我已经阅读了很多描述它的数学资料和动画代码，但我正在努力实现我的实现。

我目前的代码如下【应该可以单独在python notebook中运行】

import matplotlib.pyplot as plt 
import scipy
import cmath 
import math
import numpy as np
from matplotlib import collections  as mc
import pylab as pl


midpoints = [[305.0, 244.5], [293.5, 237.0], [274.5, 367.5], [270.5, 373.5], [229.5, 391.0], [216.0, 396.0], [302.0, 269.0], [295.0, 271.0], [60.5, 146.5], [54.0, 153.0], [52.0, 167.0], [54.0, 153.0], [52.0, 167.0], [45.0, 178.0], [75.0, 76.5], [68.5, 98.5], [75.0, 76.5], [97.0, 58.5], [283.5, 357.5], [274.5, 367.5], [309.0, 255.0], [305.0, 244.5], [309.0, 255.0], [302.0, 269.0], [299.5, 291.5], [300.0, 297.0], [300.0, 297.0], [295.0, 309.5], [62.5, 105.0], [61.0, 118.5], [62.5, 105.0], [68.5, 98.5], [58.0, 139.0], [60.5, 146.5], [241.0, 111.5], [252.0, 124.5], [256.0, 132.5], [252.0, 124.5], [283.0, 356.0], [283.5, 357.5], [300.0, 290.5], [299.5, 291.5], [300.0, 290.5], [296.0, 280.5], [296.0, 280.5], [295.0, 271.0], [158.0, 387.0], [177.0, 396.5], [197.5, 402.0], [192.5, 403.0], [189.5, 400.5], [192.5, 403.0], [197.5, 402.0], [202.5, 401.0], [214.0, 395.5], [216.0, 396.0], [202.5, 401.0], [214.0, 395.5], [233.5, 375.0], [229.5, 391.0], [233.5, 375.0], [249.0, 372.5], [282.5, 340.0], [284.5, 328.0], [284.5, 328.0], [295.0, 309.5], [45.0, 178.0], [49.5, 189.0], [57.0, 198.0], [49.5, 189.0], [238.5, 108.5], [241.0, 111.5], [162.0, 57.5], [170.0, 59.0], [239.5, 204.5], [239.0, 200.0], [293.5, 237.0], [291.0, 227.5], [265.0, 229.5], [291.0, 227.5], [239.0, 189.5], [245.5, 178.0], [239.0, 189.5], [241.0, 193.5], [241.0, 193.5], [239.0, 200.0], [55.0, 119.0], [61.0, 118.5], [53.5, 134.0], [58.0, 139.0], [50.0, 129.0], [55.0, 119.0], [50.0, 129.0], [53.5, 134.0], [107.0, 46.0], [119.5, 50.5], [97.0, 54.0], [97.0, 58.5], [107.0, 46.0], [97.0, 54.0], [150.5, 377.0], [158.0, 387.0], [257.5, 367.5], [270.5, 373.5], [249.0, 372.5], [257.5, 367.5], [280.0, 349.5], [282.5, 340.0], [280.0, 349.5], [283.0, 356.0], [239.5, 90.0], [238.5, 98.0], [238.5, 108.5], [238.5, 98.0], [130.0, 49.0], [119.5, 50.5], [189.0, 65.0], [191.0, 64.5], [189.0, 65.0], [177.0, 62.5], [170.0, 59.0], [177.0, 62.5], [256.0, 132.5], [257.5, 139.5], [128.0, 361.5], [127.5, 360.0], [136.5, 382.5], [131.5, 378.5], [126.5, 370.0], [131.5, 378.5], [128.0, 361.5], [126.5, 370.0], [105.5, 343.5], [101.0, 324.5], [105.5, 343.5], [121.5, 347.5], [126.0, 353.0], [127.5, 360.0], [121.5, 347.5], [126.0, 353.0], [191.0, 64.5], [198.5, 72.0], [237.5, 83.5], [239.5, 90.0], [145.5, 49.0], [138.5, 49.0], [159.0, 57.0], [162.0, 57.5], [145.5, 49.0], [159.0, 57.0], [265.0, 229.5], [254.5, 220.0], [253.0, 216.5], [254.5, 220.0], [253.0, 216.5], [248.0, 208.5], [248.0, 208.5], [239.5, 204.5], [245.0, 173.5], [245.5, 178.0], [250.0, 158.0], [245.0, 173.5], [257.5, 139.5], [250.0, 158.0], [177.0, 396.5], [181.0, 395.5], [181.0, 395.5], [189.5, 400.5], [147.0, 377.0], [150.5, 377.0], [140.5, 381.5], [147.0, 377.0], [140.5, 381.5], [136.5, 382.5], [92.5, 313.5], [101.0, 324.5], [99.5, 290.0], [92.5, 313.5], [98.0, 271.0], [99.5, 290.0], [134.5, 47.5], [130.0, 49.0], [134.5, 47.5], [138.5, 49.0], [73.0, 222.5], [71.0, 214.5], [107.5, 246.0], [110.5, 248.0], [104.0, 266.5], [98.0, 271.0], [69.0, 209.5], [71.0, 214.5], [226.5, 87.0], [237.5, 83.5], [226.5, 87.0], [205.5, 94.5], [205.5, 94.5], [205.5, 77.5], [198.5, 72.0], [205.5, 77.5], [96.5, 244.5], [107.5, 246.0], [109.0, 265.5], [104.0, 266.5], [108.0, 263.5], [110.5, 248.0], [109.0, 265.5], [108.0, 263.5], [65.0, 205.5], [69.0, 209.5], [57.0, 198.0], [56.0, 201.5], [65.0, 205.5], [56.0, 201.5], [90.0, 240.0], [96.5, 244.5], [83.0, 224.0], [73.0, 222.5], [90.5, 231.5], [83.0, 224.0], [90.0, 240.0], [90.5, 231.5]]

x_list = [ p[0] for p in midpoints ]
y_list = [ p[1] for p in midpoints ]
complexmdpts = [ [p[0]+1j*p[1]] for p in midpoints ]

plt.scatter(x_list, y_list, s=50, marker="x", color='y') 

coefs = np.fft.fftshift(scipy.fft.fft(complexmdpts))
n = len(coefs)
print("coeffs[{}]:\n{}".format(n, coefs[:5]))

#todo: sort coeffs?

# function in terms of t to trace out curve
def f(t):
    ftx=0
    fty=0
    for i in range(-int(n/2), int(n/2)+1):
        ftx+=(coefs[i]*cmath.exp(1j*2*math.pi*i*t/n)).real.tolist()[0]
        fty+=(coefs[i]*cmath.exp(1j*2*math.pi*i*t/n)).imag.tolist()[0]
    return [ftx/n, fty/n]
    

lines = [] # store computed lines segments to approximate function

pft = f(0) # compute first point

t_list = np.linspace(0, 2*math.pi, n) # compute list of dts to use when drawing

for t in t_list: 
    cft = f(t)
#     print("f({}): {} \n".format(t, cft))
    lines.append([cft,pft])
    pft = cft

#draw f(t) approximation
lc = mc.LineCollection(lines)
fig, ax = pl.subplots()
ax.add_collection(lc)
ax.autoscale()
ax.margins(0.1)

这是我的输出：

我想概述的要点

我的函数逼近

我不相信正确使用快速傅立叶变换。从我读到的 fft 是我需要的，我移动它是因为 scipy fft 返回移动的数组，我认为我的代码的其余部分是正确的，假设系数是正确的，这就是我怀疑的原因系数。

变换和我缺少的系数之间是否有一个步骤？或者我的函数评估是否给定系数不正确？还是我错过了其他东西？

Answer 1

实际上有一些问题需要解决才能获得预期的大纲。让我们逐一讨论这些问题。

傅立叶系数计算

由于您正在计算 2D 数组的 FFT complexmdpts（每个元素都是 1 个复数值的数组），因此默认行为scipy.fft是沿最后一个轴计算 FFT。在这种情况下，这意味着您实际上正在计算n长度为 1 的 FFT，并且鉴于长度为 1 的 FFT 是恒等式，整个计算将返回原始数组。一种解决方案是axis=0明确指定：

scipy.fft(complexmdpts,axis=0) 

鉴于基本上只有一维，您也可以将二维数组展平为一维数组，但这会导致对基于此二维结构的其余代码进行更多更改。

系数偏移

在尝试解释 FFT 系数时，通常会出现混淆。较高指数中的系数对应于发生在奈奎斯特频率之上的环绕之后的负频率。为了在负频率出现在其直观位置的图上可视化这些系数，将较高索引中fftshift的系数与较低索引中的系数交换。这样，对应于负频率的系数出现在对应于正频率的系数之前。

在您的特定情况下，将f(t)负频率（无论何时i为负， 的参数也是如此cmath.exp）与负索引处的系数相关联的计算。由于 python 数组索引方便地使用从数组末尾倒数的元素作为负索引，因此在使用负索引进行索引时，您正确地使用了与负频率相对应的系数。所有这一切都说明您不需要用 交换系数fftshift，并且可以直接通过以下方式获得系数：

coefs = scipy.fft(complexmdpts,axis=0)

采样域

您已指定t范围从 0 到2*math.pi. 鉴于您f(t)的域实现实际上范围从 0 到n（例如，当i==1参数 tocmath.exp从 0 到1j*2*math.piast从 0 到 时n）。要使曲线跨越整个域，您应该t相应地更新您的值：

t_list = np.linspace(0, n, n)

积分排序

最后，通过使用散点图绘制原始点系列，您隐藏了线图的外观：

试图获得一个平滑的函数来插值这个结果，结果是一个相似的线条从轮廓的一侧跳到另一侧。如果你想捕捉这些跳跃，那么我想你已经完成了。但我怀疑您可能想要捕捉该区域的轮廓。在这种情况下，您必须对原始点进行排序，以便它们遵循大纲。一种方法是迭代地将最近点附加到先前选择的点。这将给出更接近轮廓的东西：

出于演示目的（即没有声称这是最好的方法），您可以使用以下内容进行排序：

def sort_points(points):
  # pick a point
  reference_point = points[0]
  sorted = [reference_point]
  remaining_points = range(1,len(points))
  for i in range(1,len(points)):

    # find the closest point to reference_point, 
    mindiff = np.sum(np.square(np.array(points[remaining_points[0]])-reference_point))
    idx = 0
    # loop over all the other remaining points
    for j in range(1,len(remaining_points)):
      diff = np.sum(np.square(np.array(points[remaining_points[j]])-reference_point))
      if diff < mindiff:
        mindiff = diff
        idx = j        
    # found the closest: update the selected point, and add it to the list of sorted points
    reference_point = points[remaining_points[idx]]
    sorted.append(reference_point )
    remaining_points = np.delete(remaining_points, idx)
  return sorted    

作为参考，这里是完整的代码：

import matplotlib.pyplot as plt 
import scipy
import cmath 
import math
import numpy as np
from matplotlib import collections  as mc
import pylab as pl

def sort_points(points):
  # pick a point
  reference_point = points[0]
  sorted = [reference_point]
  remaining_points = range(1,len(points))
  for i in range(1,len(points)):

    # find the closest point to reference_point, 
    mindiff = np.sum(np.square(np.array(points[remaining_points[0]])-reference_point))
    idx = 0
    # loop over all the other remaining points
    for j in range(1,len(remaining_points)):
      diff = np.sum(np.square(np.array(points[remaining_points[j]])-reference_point))
      if diff < mindiff:
        mindiff = diff
        idx = j        
    # found the closest: update the selected point, and add it to the list of sorted points
    reference_point = points[remaining_points[idx]]
    sorted.append(reference_point )
    remaining_points = np.delete(remaining_points, idx)
  return sorted    


midpoints = [[305.0, 244.5], [293.5, 237.0], [274.5, 367.5], [270.5, 373.5], [229.5, 391.0], [216.0, 396.0], [302.0, 269.0], [295.0, 271.0], [60.5, 146.5], [54.0, 153.0], [52.0, 167.0], [54.0, 153.0], [52.0, 167.0], [45.0, 178.0], [75.0, 76.5], [68.5, 98.5], [75.0, 76.5], [97.0, 58.5], [283.5, 357.5], [274.5, 367.5], [309.0, 255.0], [305.0, 244.5], [309.0, 255.0], [302.0, 269.0], [299.5, 291.5], [300.0, 297.0], [300.0, 297.0], [295.0, 309.5], [62.5, 105.0], [61.0, 118.5], [62.5, 105.0], [68.5, 98.5], [58.0, 139.0], [60.5, 146.5], [241.0, 111.5], [252.0, 124.5], [256.0, 132.5], [252.0, 124.5], [283.0, 356.0], [283.5, 357.5], [300.0, 290.5], [299.5, 291.5], [300.0, 290.5], [296.0, 280.5], [296.0, 280.5], [295.0, 271.0], [158.0, 387.0], [177.0, 396.5], [197.5, 402.0], [192.5, 403.0], [189.5, 400.5], [192.5, 403.0], [197.5, 402.0], [202.5, 401.0], [214.0, 395.5], [216.0, 396.0], [202.5, 401.0], [214.0, 395.5], [233.5, 375.0], [229.5, 391.0], [233.5, 375.0], [249.0, 372.5], [282.5, 340.0], [284.5, 328.0], [284.5, 328.0], [295.0, 309.5], [45.0, 178.0], [49.5, 189.0], [57.0, 198.0], [49.5, 189.0], [238.5, 108.5], [241.0, 111.5], [162.0, 57.5], [170.0, 59.0], [239.5, 204.5], [239.0, 200.0], [293.5, 237.0], [291.0, 227.5], [265.0, 229.5], [291.0, 227.5], [239.0, 189.5], [245.5, 178.0], [239.0, 189.5], [241.0, 193.5], [241.0, 193.5], [239.0, 200.0], [55.0, 119.0], [61.0, 118.5], [53.5, 134.0], [58.0, 139.0], [50.0, 129.0], [55.0, 119.0], [50.0, 129.0], [53.5, 134.0], [107.0, 46.0], [119.5, 50.5], [97.0, 54.0], [97.0, 58.5], [107.0, 46.0], [97.0, 54.0], [150.5, 377.0], [158.0, 387.0], [257.5, 367.5], [270.5, 373.5], [249.0, 372.5], [257.5, 367.5], [280.0, 349.5], [282.5, 340.0], [280.0, 349.5], [283.0, 356.0], [239.5, 90.0], [238.5, 98.0], [238.5, 108.5], [238.5, 98.0], [130.0, 49.0], [119.5, 50.5], [189.0, 65.0], [191.0, 64.5], [189.0, 65.0], [177.0, 62.5], [170.0, 59.0], [177.0, 62.5], [256.0, 132.5], [257.5, 139.5], [128.0, 361.5], [127.5, 360.0], [136.5, 382.5], [131.5, 378.5], [126.5, 370.0], [131.5, 378.5], [128.0, 361.5], [126.5, 370.0], [105.5, 343.5], [101.0, 324.5], [105.5, 343.5], [121.5, 347.5], [126.0, 353.0], [127.5, 360.0], [121.5, 347.5], [126.0, 353.0], [191.0, 64.5], [198.5, 72.0], [237.5, 83.5], [239.5, 90.0], [145.5, 49.0], [138.5, 49.0], [159.0, 57.0], [162.0, 57.5], [145.5, 49.0], [159.0, 57.0], [265.0, 229.5], [254.5, 220.0], [253.0, 216.5], [254.5, 220.0], [253.0, 216.5], [248.0, 208.5], [248.0, 208.5], [239.5, 204.5], [245.0, 173.5], [245.5, 178.0], [250.0, 158.0], [245.0, 173.5], [257.5, 139.5], [250.0, 158.0], [177.0, 396.5], [181.0, 395.5], [181.0, 395.5], [189.5, 400.5], [147.0, 377.0], [150.5, 377.0], [140.5, 381.5], [147.0, 377.0], [140.5, 381.5], [136.5, 382.5], [92.5, 313.5], [101.0, 324.5], [99.5, 290.0], [92.5, 313.5], [98.0, 271.0], [99.5, 290.0], [134.5, 47.5], [130.0, 49.0], [134.5, 47.5], [138.5, 49.0], [73.0, 222.5], [71.0, 214.5], [107.5, 246.0], [110.5, 248.0], [104.0, 266.5], [98.0, 271.0], [69.0, 209.5], [71.0, 214.5], [226.5, 87.0], [237.5, 83.5], [226.5, 87.0], [205.5, 94.5], [205.5, 94.5], [205.5, 77.5], [198.5, 72.0], [205.5, 77.5], [96.5, 244.5], [107.5, 246.0], [109.0, 265.5], [104.0, 266.5], [108.0, 263.5], [110.5, 248.0], [109.0, 265.5], [108.0, 263.5], [65.0, 205.5], [69.0, 209.5], [57.0, 198.0], [56.0, 201.5], [65.0, 205.5], [56.0, 201.5], [90.0, 240.0], [96.5, 244.5], [83.0, 224.0], [73.0, 222.5], [90.5, 231.5], [83.0, 224.0], [90.0, 240.0], [90.5, 231.5]]
midpoints = sort_points(midpoints)

x_list = [ p[0] for p in midpoints ]
y_list = [ p[1] for p in midpoints ]
complexmdpts = [ [p[0]+1j*p[1]] for p in midpoints ]

plt.scatter(x_list, y_list, s=50, marker="x", color='y') 

coefs = scipy.fft(complexmdpts,axis=0)
n = len(coefs)
print("coeffs[{}]:\n{}".format(n, coefs[:5]))

# function in terms of t to trace out curve
m = n
def f(t):
    ftx=0
    fty=0
    for i in range(-int(m/2), int(m/2)+1):
        ftx+=(coefs[i]*cmath.exp(1j*2*math.pi*i*t/n)).real.tolist()[0]
        fty+=(coefs[i]*cmath.exp(1j*2*math.pi*i*t/n)).imag.tolist()[0]
    return [ftx/n, fty/n]
    

lines = [] # store computed lines segments to approximate function

pft = f(0) # compute first point

t_list = np.linspace(0, n, n) # compute list of dts to use when drawing

for t in t_list: 
    cft = f(t)
    lines.append([cft,pft])
    pft = cft

#draw f(t) approximation
lc = mc.LineCollection(lines)
fig, ax = pl.subplots()
ax.add_collection(lc)
ax.autoscale()
ax.margins(0.1)