Python 如何在Numpy中计算傅立叶级数？
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/4258106/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
How to calculate a Fourier series in Numpy?
提问by Mermoz
I have a periodic function of period T and would like to know how to obtain the list of the Fourier coefficients. I tried using fftmodule from numpy but it seems more dedicated to Fourier transforms than series. Maybe it a lack of mathematical knowledge, but I can't see how to calculate the Fourier coefficients from fft.
我有一个周期为 T 的周期函数，想知道如何获得傅立叶系数列表。我尝试使用numpy 的fft模块，但它似乎更专注于傅立叶变换而不是系列。也许是缺乏数学知识，但我看不出如何从 fft 计算傅立叶系数。
Help and/or examples appreciated.
帮助和/或示例表示赞赏。
采纳答案by Mermoz
In the end, the most simple thing (calculating the coefficient with a riemann sum) was the most portable/efficient/robust way to solve my problem:
最后，最简单的事情（用黎曼和计算系数）是解决我的问题的最便携/最有效/最稳健的方法：
def cn(n):
c = y*np.exp(-1j*2*n*np.pi*time/period)
return c.sum()/c.size
def f(x, Nh):
f = np.array([2*cn(i)*np.exp(1j*2*i*np.pi*x/period) for i in range(1,Nh+1)])
return f.sum()
y2 = np.array([f(t,50).real for t in time])
plot(time, y)
plot(time, y2)
gives me:
给我：
回答by mtrw
Do you have a list of discrete samples of your function, or is your function itself discrete? If so, the Discrete Fourier Transform, calculated using an FFT algorithm, provides the Fourier coefficients directly (see here).
您有函数的离散样本列表，还是函数本身是离散的？如果是这样，则使用 FFT 算法计算的离散傅立叶变换会直接提供傅立叶系数（请参阅此处）。
On the other hand, if you have an analytic expression for the function, you probably need a symbolic math solver of some kind.
另一方面，如果您有函数的解析表达式，您可能需要某种符号数学求解器。
回答by dr jimbob
Numpy isn't the right tool really to calculate fourier series components, as your data has to be discretely sampled. You really want to use something like Mathematica or should be using fourier transforms.
Numpy 并不是真正计算傅立叶级数分量的正确工具，因为您的数据必须是离散采样的。你真的想使用像 Mathematica 这样的东西，或者应该使用傅立叶变换。
To roughly do it, let's look at something simple a triangle wave of period 2pi, where we can easily calculate the Fourier coefficients (c_n = -i ((-1)^(n+1))/n for n>0; e.g., c_n = { -i, i/2, -i/3, i/4, -i/5, i/6, ... } for n=1,2,3,4,5,6 (using Sum( c_n exp(i 2 pi n x) ) as Fourier series).
为了粗略地做到这一点，让我们看一些简单的周期为 2pi 的三角波，我们可以很容易地计算出傅立叶系数 (c_n = -i ((-1)^(n+1))/n for n>0；例如, c_n = { -i, i/2, -i/3, i/4, -i/5, i/6, ... } for n=1,2,3,4,5,6（使用 Sum ( c_n exp(i 2 pi nx) ) 作为傅立叶级数)。
import numpy
x = numpy.arange(0,2*numpy.pi, numpy.pi/1000)
y = (x+numpy.pi/2) % numpy.pi - numpy.pi/2
fourier_trans = numpy.fft.rfft(y)/1000
If you look at the first several Fourier components:
如果您查看前几个傅立叶分量：
array([ -3.14159265e-03 +0.00000000e+00j,
2.54994550e-16 -1.49956612e-16j,
3.14159265e-03 -9.99996710e-01j,
1.28143395e-16 +2.05163971e-16j,
-3.14159265e-03 +4.99993420e-01j,
5.28320925e-17 -2.74568926e-17j,
3.14159265e-03 -3.33323464e-01j,
7.73558750e-17 -3.41761974e-16j,
-3.14159265e-03 +2.49986840e-01j,
1.73758496e-16 +1.55882418e-17j,
3.14159265e-03 -1.99983550e-01j,
-1.74044469e-16 -1.22437710e-17j,
-3.14159265e-03 +1.66646927e-01j,
-1.02291982e-16 -2.05092972e-16j,
3.14159265e-03 -1.42834113e-01j,
1.96729377e-17 +5.35550532e-17j,
-3.14159265e-03 +1.24973680e-01j,
-7.50516717e-17 +3.33475329e-17j,
3.14159265e-03 -1.11081501e-01j,
-1.27900121e-16 -3.32193126e-17j,
-3.14159265e-03 +9.99670992e-02j,
First neglect the components that are near 0 due to floating point accuracy (~1e-16, as being zero). The more difficult part is to see that the 3.14159 numbers (that arose before we divide by the period of a 1000) should also be recognized as zero, as the function is periodic). So if we neglect those two factors we get:
首先忽略由于浮点精度而接近 0 的分量（~1e-16，为零）。更困难的部分是看到 3.14159 数字（在我们除以 1000 的周期之前出现的）也应该被识别为零，因为函数是周期性的）。因此，如果我们忽略这两个因素，我们会得到：
fourier_trans = [0,0,-i,0,i/2,0,-i/3,0,i/4,0,-i/5,0,-i/6, ...
and you can see the fourier series numbers come up as every other number (I haven't investigated; but I believe the components correspond to [c0, c-1, c1, c-2, c2, ... ]). I'm using conventions according to wiki: http://en.wikipedia.org/wiki/Fourier_series.
并且您可以看到傅立叶级数与其他数字一样出现（我还没有研究过；但我相信这些分量对应于 [c0, c-1, c1, c-2, c2, ... ]）。我根据 wiki 使用约定：http: //en.wikipedia.org/wiki/Fourier_series。
Again, I'd suggest using mathematica or a computer algebra system capable of integrating and dealing with continuous functions.
同样，我建议使用 mathematica 或能够集成和处理连续函数的计算机代数系统。
回答by DaveP
As other answers have mentioned, it seems that what you are looking for is a symbolic computing package, so numpy isn't suitable. If you wish to use a free python-based solution, then either sympyor sageshould meet your needs.
正如其他答案所提到的，您正在寻找的是一个符号计算包，因此 numpy 不合适。如果您希望使用免费的基于 Python 的解决方案，那么sympy或sage应该可以满足您的需求。
回答by gg349
This is an old question, but since I had to code this, I am posting here the solution that uses the numpy.fftmodule, that is likely faster than other hand-crafted solutions.
这是一个老问题，但由于我必须对此进行编码，因此我在此处发布了使用该numpy.fft模块的解决方案，这可能比其他手工制作的解决方案更快。
The DFTis the right tool for the job of calculating up to numerical precision the coefficients of the Fourier series of a function, defined as an analytic expression of the argument or as a numerical interpolating function over some discrete points.
该DFT是用于计算达到数值精度的傅立叶级数的函数，定义为参数的解析表达式或作为以上一些离散的点的数值进行内插函数的系数的工作的工具。
This is the implementation, which allows to calculate the real-valued coefficients of the Fourier series, or the complex valued coefficients, by passing an appropriate return_complex:
这是实现，它允许通过传递适当的 来计算傅立叶级数的实值系数或复值系数return_complex：
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
This is an example of usage:
这是一个使用示例：
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])
And the plot of the resulting a0,a1,...,a10,b1,b2,...,b10coefficients:
以及所得a0,a1,...,a10,b1,b2,...,b10系数的图：
This is an optional test for the function, for both modes of operation. You should run this after the example, or define a periodic function fand a period Tbefore running the code.
这是功能的可选测试，适用于两种操作模式。您应该在示例之后运行它，或者在运行代码之前定义一个周期函数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)
