这是一个老问题，但是由于我必须编写代码，所以我在这里发布了使用numpy.fft模块的解决方案，这可能比其他手工编制的解决方案更快。

DFT是计算函数Fourier级数系数(定义为参数的解析表达式或某些离散点上的数值插值函数)达到数值精度的正确工具。

这是一个实现，它允许通过传递适当的return_complex来计算Fourier级数的实值系数或复数系数：def fourier_series_coeff_numpy(f, T, N, return_complex=False):

"""Calculates the first 2*N+1 Fourier series coeff. of a periodic function.

Given a periodic, function f(t) with period T, this function returns the

coefficients a0, {a1,a2,...},{b1,b2,...} such that:

f(t) ~= a0/2+ sum_{k=1}^{N} ( a_k*cos(2*pi*k*t/T) + b_k*sin(2*pi*k*t/T) )

If return_complex is set to True, it returns instead the coefficients

{c0,c1,c2,...}

such that:

f(t) ~= sum_{k=-N}^{N} c_k * exp(i*2*pi*k*t/T)

where we define c_{-n} = complex_conjugate(c_{n})

Refer to wikipedia for the relation between the real-valued and complex

valued coeffs at http://en.wikipedia.org/wiki/Fourier_series.

Parameters

----------

f : the periodic function, a callable like f(t)

T : the period of the function f, so that f(0)==f(T)

N_max : the function will return the first N_max + 1 Fourier coeff.

Returns

-------

if return_complex == False, the function returns:

a0 : float

a,b : numpy float arrays describing respectively the cosine and sine coeff.

if return_complex == True, the function returns:

c : numpy 1-dimensional complex-valued array of size N+1

"""

# From Shanon theoreom we must use a sampling freq. larger than the maximum

# frequency you want to catch in the signal.

f_sample = 2 * N

# we also need to use an integer sampling frequency, or the

# points will not be equispaced between 0 and 1. We then add +2 to f_sample

t, dt = np.linspace(0, T, f_sample + 2, endpoint=False, retstep=True)

y = np.fft.rfft(f(t)) / t.size

if return_complex:

return y

else:

y *= 2

return y[0].real, y[1:-1].real, -y[1:-1].imag

这是一个用法示例：from numpy import ones_like, cos, pi, sin, allclose

T = 1.5 # any real number

def f(t):

"""example of periodic function in [0,T]"""

n1, n2, n3 = 1., 4., 7. # in Hz, or nondimensional for the matter.

a0, a1, b4, a7 = 4., 2., -1., -3

return a0 / 2 * ones_like(t) + a1 * cos(2 * pi * n1 * t / T) + b4 * sin(

2 * pi * n2 * t / T) + a7 * cos(2 * pi * n3 * t / T)

N_chosen = 10

a0, a, b = fourier_series_coeff_numpy(f, T, N_chosen)

# we have as expected that

assert allclose(a0, 4)

assert allclose(a, [2, 0, 0, 0, 0, 0, -3, 0, 0, 0])

assert allclose(b, [0, 0, 0, -1, 0, 0, 0, 0, 0, 0])

以及由此产生的a0,a1,...,a10,b1,b2,...,b10系数的图：

对于两种操作模式，这是功能的可选测试。您应该在示例之后运行此函数，或者在运行代码之前定义周期函数f和周期T。# #### test that it works with real coefficients:

from numpy import linspace, allclose, cos, sin, ones_like, exp, pi, \

complex64, zeros

def series_real_coeff(a0, a, b, t, T):

"""calculates the Fourier series with period T at times t,

from the real coeff. a0,a,b"""

tmp = ones_like(t) * a0 / 2.

for k, (ak, bk) in enumerate(zip(a, b)):

tmp += ak * cos(2 * pi * (k + 1) * t / T) + bk * sin(

2 * pi * (k + 1) * t / T)

return tmp

t = linspace(0, T, 100)

f_values = f(t)

a0, a, b = fourier_series_coeff_numpy(f, T, 52)

# construct the series:

f_series_values = series_real_coeff(a0, a, b, t, T)

# check that the series and the original function match to numerical precision:

assert allclose(f_series_values, f_values, atol=1e-6)

# #### test similarly that it works with complex coefficients:

def series_complex_coeff(c, t, T):

"""calculates the Fourier series with period T at times t,

from the complex coeff. c"""

tmp = zeros((t.size), dtype=complex64)

for k, ck in enumerate(c):

# sum from 0 to +N

tmp += ck * exp(2j * pi * k * t / T)

# sum from -N to -1

if k != 0:

tmp += ck.conjugate() * exp(-2j * pi * k * t / T)

return tmp.real

f_values = f(t)

c = fourier_series_coeff_numpy(f, T, 7, return_complex=True)

f_series_values = series_complex_coeff(c, t, T)

assert allclose(f_series_values, f_values, atol=1e-6)