如何在Python中用FFT计算傅里叶系数并匹配基频

我有一条XYZ坐标系下的曲线，希望将其按指定形式进行傅里叶展开。我想通过Python中的FFT计算展开式中的$x_{c,m}$、$x_{s,m}$等系数，但在匹配基频与$\theta$以获取正确系数时遇到困难。我需要的是该特定形式的展开，而非其他类型。

当前使用的代码如下：

import numpy as np
from scipy import interpolate
from scipy.fft import rfft, rfftfreq
from math import pi
import matplotlib.pyplot as plt

def fourier_series(a0,a,b,order,x):
    return a0 + np.sum([a[n]*np.cos(n*x) + b[n]*np.sin(n*x) for n in range(order)])

order = 25
theta = np.linspace(0,2*pi,100)
curve = [[np.cos(x), 2*np.sin(x), np.cos(2*x) + np.sin(2*x)] for x in theta] #this could be any closed curve.

xArr, yArr, zArr = np.transpose(curve)

# Calculate the fft
freq = rfftfreq(len(xArr))
freq_series_cos = [n/len(xArr) for n in range(1,order+1)]
freq_series_sin = [n/len(xArr) for n in range(order+1)]
curvesFourier = []

#interpolate the fft and calculate the right frequencies to match CurveXYZFourier
for x in [xArr, yArr, zArr]:
    xf = rfft(x)/len(x)
    fft_0 = pi*xf[0].real
    fft_cos = pi*xf.real/2  #find the cosine coefficients
    fft_sin = -pi*xf.imag/2  #find the sine coefficients
    
    fft_cos_interpolated = interpolate.CubicSpline(freq,fft_cos) #interpolate the coefficients to pick the right frequencies 
    fft_sin_interpolated = interpolate.CubicSpline(freq,fft_sin) #interpolate the coefficients to pick the right frequencies
    
    b  = fft_sin_interpolated(freq_series_sin) 
    a0 = fft_0 
    a  = fft_cos_interpolated(freq_series_cos)
    
    x_fourier = [fourier_series(a0,a,b,order, i) for i in theta]
    curvesFourier.append(x_fourier)


ax = plt.axes(projection='3d')

ax.plot3D(xArr, yArr, zArr, "k--")
ax.scatter3D(curvesFourier[0], curvesFourier[1], curvesFourier[2])
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

plt.show()

目前我通过插值频率来处理，但我认为应该有无需插值、直接匹配目标基频的方法。我也曾尝试适配Stack Overflow上gg349的解决方案，但未成功。

核心问题在于对FFT频率点的映射理解有误——不需要插值，FFT输出的频率点正好对应傅里叶级数的阶数，直接索引即可。

以下是修正后的代码，去掉冗余插值步骤，直接匹配基频与系数：

import numpy as np
from scipy.fft import rfft, rfftfreq
from math import pi
import matplotlib.pyplot as plt

def fourier_series(a0, a, b, order, x):
    # 修正索引：a和b对应1~order阶，循环从1开始
    series = a0
    for n in range(1, order+1):
        series += a[n-1] * np.cos(n*x) + b[n-1] * np.sin(n*x)
    return series

order = 25
# 去掉端点避免重复采样，保证FFT频率映射准确
theta = np.linspace(0, 2*pi, 100, endpoint=False)
curve = [[np.cos(x), 2*np.sin(x), np.cos(2*x) + np.sin(2*x)] for x in theta]

xArr, yArr, zArr = np.transpose(curve)
curvesFourier = []

for coord in [xArr, yArr, zArr]:
    xf = rfft(coord)
    N = len(coord)
    
    # 直接从FFT结果提取对应阶数的系数
    a0 = (xf[0].real) / N * 2 * pi  # 直流分量匹配展开式缩放
    a = (xf[1:order+1].real) / N * pi  # 余弦系数对应FFT实部
    b = -(xf[1:order+1].imag) / N * pi  # 正弦系数对应FFT虚部的相反数
    
    # 生成拟合曲线
    x_fourier = [fourier_series(a0, a, b, order, t) for t in theta]
    curvesFourier.append(x_fourier)

# 绘图对比
ax = plt.axes(projection='3d')
ax.plot3D(xArr, yArr, zArr, "k--", label="Original Curve")
ax.plot3D(curvesFourier[0], curvesFourier[1], curvesFourier[2], "r-", label="Fourier Fit")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.legend()
plt.show()

关键修正说明：

• 采样优化：用endpoint=False生成θ，避免0和2π重复采样，确保FFT频率映射精准。
• 直接索引系数：FFT输出的第m个点（从1开始）对应傅里叶级数的m阶分量，直接取xf[1:order+1]即可，无需插值。
• 系数缩放匹配：根据傅里叶级数与FFT的数学关系调整缩放因子，完全匹配你需要的展开形式。
• 级数计算修正：修复原函数的循环索引错误，确保系数与阶数一一对应。

处理后拟合曲线会与原曲线几乎完全重合，且无需任何插值操作，直接通过索引完成基频与系数的匹配。

内容的提问来源于stack exchange，提问作者Miguel Madeira