当前位置： 文江博客话题详情

如何在 Numpy 中计算傅立叶级数？

发布于 2024-10-04 17:25:41 字数 233 浏览 8 评论 0原文

我有一个周期为 T 的周期函数，想知道如何获取傅立叶系数列表。我尝试使用 numpy 中的 fft 模块，但它似乎更专注于傅里叶变换而不是级数。 也许是缺乏数学知识，但我不知道如何从 fft 计算傅里叶系数。

感谢帮助和/或示例。

原文

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 fft module 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.

Help and/or examples appreciated.

如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论， 你可以免费 注册 一个本站的账号。

评论（5）

谁的新欢旧爱 2024-10-11 17:25:41

最后，最简单的事情（用黎曼和计算系数）是解决我的问题的最便携/有效/稳健的方法：

import numpy as np
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)

给我：
替代文本

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:

import numpy as np
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:
alt text

回复 原文
罪歌 2024-10-11 17:25:41

这是一个老问题，但由于我必须对此进行编码，因此我在这里发布使用 numpy.fft 模块的解决方案，这可能比其他手工制作的解决方案更快。

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

这是一个实现，它允许通过传递适当的 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

这是一个用法示例：

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)

This is an old question, but since I had to code this, I am posting here the solution that uses the numpy.fft module, that is likely faster than other hand-crafted solutions.

The DFT is 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.

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:

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,...,b10 coefficients:
enter image description here

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 f and a period T before running the code.

# #### 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)
回复 原文
水波映月 2024-10-11 17:25:41

Numpy 并不是真正计算傅里叶级数分量的正确工具，因为您的数据必须进行离散采样。您确实想使用 Mathematica 之类的东西，或者应该使用傅里叶变换。

为了粗略地做到这一点，让我们看一些简单的周期为 2pi 的三角波，其中我们可以轻松计算傅里叶系数（c_n = -i ((-1)^(n+1))/n for n>0; 例如, c_n = { -i, i/2, -i/3, i/4, -i/5, i/6, ... } 对于 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

如果您查看前几个傅立叶分量：

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,

首先忽略由于浮点精度而接近 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, ...

你可以看到傅里叶级数与其他数字一样出现（我没有调查；但我相信这些分量对应于 [c0, c-1, c1, c-2 ，c2，...]）。我根据 wiki 使用约定： http://en.wikipedia.org/wiki/Fourier_series。

再次，我建议使用能够积分和处理连续函数的mathematica 或计算机代数系统。

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.

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

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:

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.

Again, I'd suggest using mathematica or a computer algebra system capable of integrating and dealing with continuous functions.

回复 原文
怪异←思 2024-10-11 17:25:41

正如其他答案所提到的，您似乎正在寻找的是符号计算包，因此 numpy 不适合。如果您希望使用基于 Python 的免费解决方案，则可以使用 sympy 或 sage 应该可以满足您的需求。

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 sympy or sage should meet your needs.

回复 原文
错々过的事 2024-10-11 17:25:41

您是否有函数的离散样本列表，或者您的函数本身是离散的？如果是这样，使用 FFT 算法计算的离散傅里叶变换将直接提供傅里叶系数（请参阅此处）。

另一方面，如果您有该函数的解析表达式，您可能需要某种符号数学求解器。

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

On the other hand, if you have an analytic expression for the function, you probably need a symbolic math solver of some kind.

回复 原文
~没有更多了~

关于作者

泪之魂

暂无简介

文章
评论
25 人气
发私信

相关话题

更多

热门标签

操作系统 程序设计 IT运维 Linux系统管理 JavaScript 服务器应用 solaris C/C++ PHP Shell BSD Vue.js aix Oracle Python HTML 系统管理 HTML5 CSS 前端
更多

推荐作者

峩卟喜欢

文章 0 评论 0

一向肩并

文章 0 评论 0

潜伏

文章 0 评论 0

dongyinghao

文章 0 评论 0

百变从容

文章 0 评论 0

沧笙踏歌

文章 0 评论 0

更多

友情链接

文江博客
    我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站，即表示您同意使用 Cookies 和您的相关数据。
    原文