scipy 在Python中使用fft计算傅立叶系数并正确匹配基频

shstlldc 于 2023-05-29 发布在 Python

关注(0)|答案(1)|浏览(127)

在xyz坐标系中有一条曲线。我想把这条曲线傅立叶展开为：
Curve fourier expansion.

我想计算展开系数xc，m，xs，m，…在Python中使用FFT。我在理解如何将基频与θ匹配以便fft返回正确的系数时遇到了麻烦。请注意，我不希望使用任何扩展，我特别希望使用此扩展。
以下是我的当前代码：

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()

我在插入频率。然而，我希望有一种方法可以将基频与预期频率相匹配，因此不需要这一步。我也尝试过在https://stackoverflow.com/questions/4258106/how-to-calculate-a-fourier-series-in-numpy中采用gg 349给出的解决方案，但没有成功。
谢谢你的帮助！

scipy

来源：https://stackoverflow.com/questions/76288865/calculating-fourier-coefficients-in-python-using-fft-and-correctly-matching-the

1条答案

按热度按时间

9wbgstp71#

结果是我把问题弄得太复杂了。下面是一个解决方案：

import numpy as np
from scipy import interpolate
from scipy.fft import rfft
from math import pi
import matplotlib.pyplot as plt
   
def fourier_series(a0,a,b,order,x):
    return a0 + np.sum([a[n-1]*np.cos(n*x) for n in range(1,order)]) +np.sum([b[n]*np.sin(n*x) for n in range(order)])

order = 50
theta = np.linspace(0,2*pi,100)
curve = [[1 + 0.02*np.cos(x)+ 0.03*np.cos(2*x), 1 + 0.02*np.sin(x), np.cos(2*x) + np.sin(2*x)] for x in theta] #any closed curve

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

'''
# Compute the arclength parameterization of the curve, only needed for interpolating the data
dL = np.sqrt(np.sum(np.diff(curve, axis=0)**2, axis=1))
L = np.concatenate(([0], np.cumsum(dL)))
L = L*2*pi/L[-1]

#To find orders above len(x)/2 you must interpolate the data and choose the sampling so that len(n_samples) >= 2*order
x_interpol = interpolate.CubicSpline(L,xArr)
y_interpol = interpolate.CubicSpline(L,yArr)
z_interpol = interpolate.CubicSpline(L,zArr)
n_samples = np.linspace(0,2*pi,1000)
'''

coeffs = []
curvesFourier = []
curvesFourier2 = []

for x in [xArr, yArr, zArr]:
    #xf=rfft(x(n_samples))/len(x(n_samples)) #use this instead if you wish to find any order
    xf=rfft(x)/len(x)                        #using the real array produces more accurate results but the maximum order is limited to len(x)/2
    fft_0 = xf[0].real
    fft_cos = 2*xf[1:].real  #find the cosine coefficients
    fft_sin = -2*xf.imag  #find the sine coefficients
    x_fourier = [fourier_series(fft_0,fft_cos,fft_sin,order, i) for i in theta]
    curvesFourier.append(x_fourier)
    
ax = plt.axes(projection='3d')

ax.plot3D(xArr, yArr, zArr, "k--")
ax.plot3D(curvesFourier[0], curvesFourier[1], curvesFourier[2],"g-")

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

plt.show()

赞(0）回复(0）举报 2023-05-29

我来回答

scipy 在Python中使用fft计算傅立叶系数并正确匹配基频

1条答案

相关问题

热门标签

最新问答