计算傅立叶级数来表示数据点

如何解决计算傅立叶级数来表示数据点

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

类似于这里发生的事情 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],[45.0,178.0],[75.0,76.5],[68.5,98.5],[97.0,58.5],[283.5,357.5],[309.0,255.0],[305.0,[299.5,291.5],[300.0,297.0],309.5],[62.5,105.0],[61.0,118.5],[58.0,139.0],[241.0,111.5],[252.0,124.5],[256.0,132.5],[283.0,356.0],290.5],[296.0,280.5],[158.0,387.0],[177.0,396.5],[197.5,402.0],[192.5,403.0],[189.5,400.5],[202.5,401.0],[214.0,395.5],[233.5,375.0],[249.0,372.5],[282.5,340.0],[284.5,328.0],[49.5,189.0],[57.0,198.0],[238.5,108.5],[162.0,57.5],[170.0,59.0],[239.5,204.5],[239.0,200.0],[291.0,227.5],[265.0,229.5],189.5],[245.5,193.5],[55.0,119.0],[53.5,134.0],[50.0,129.0],[107.0,46.0],[119.5,50.5],54.0],[150.5,377.0],[257.5,[280.0,349.5],90.0],98.0],[130.0,49.0],[189.0,65.0],[191.0,64.5],62.5],139.5],[128.0,361.5],[127.5,360.0],[136.5,382.5],[131.5,378.5],[126.5,370.0],[105.5,343.5],[101.0,324.5],[121.5,347.5],[126.0,353.0],[198.5,72.0],[237.5,83.5],[145.5,[138.5,[159.0,57.0],[254.5,220.0],[253.0,216.5],[248.0,208.5],[245.0,173.5],[250.0,158.0],[181.0,[147.0,[140.5,381.5],[92.5,313.5],[99.5,290.0],[98.0,[134.5,47.5],[73.0,222.5],[71.0,214.5],[107.5,246.0],[110.5,248.0],[104.0,266.5],[69.0,209.5],[226.5,87.0],[205.5,94.5],77.5],[96.5,[109.0,265.5],[108.0,263.5],[65.0,205.5],[56.0,201.5],[90.0,240.0],[83.0,224.0],[90.5,231.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 返回移动的数组，我认为我的其余代码是正确的，假设系数是正确的，这就是为什么我怀疑系数。

变换和我遗漏的系数之间是否有一个步骤？或者我的函数评估是否给定系数不正确？还是我错过了什么？

解决方法

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

傅立叶系数计算

由于您要计算二维数组 complexmdpts 的 FFT（每个元素都是 1 个复数值的数组），scipy.fft 的默认行为是沿最后一个轴计算 FFT。在这种情况下，这意味着您实际上正在计算长度为 1 的 n 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 时 cmath.exp 的参数从 0 到 1j*2*math.pi 作为 { {1}} 从 0 到 t）。要使曲线跨越整个域，您应该相应地更新 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