Comment calculer une série de Fourier en Numpy?

Question

J'ai une fonction périodique de période T et je voudrais savoir comment obtenir la liste des coefficients de Fourier. J'ai essayé d'utiliser le module fft de numpy mais il semble plus dédié aux transformées de Fourier que les séries. Peut-être un manque de connaissances mathématiques, mais je ne vois pas comment calculer les coefficients de Fourier de fft.

Aide et/ou exemples appréciés.

Answer 1

En fin de compte, la chose la plus simple (le calcul du coefficient avec une somme de Riemann) était le moyen le plus portable/efficace/robuste pour résoudre mon problème:

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) 

me donne:

Answer 2

Avez-vous une liste d'échantillons discrets de votre fonction, ou votre fonction elle-même est-elle discrète? Si c'est le cas, la transformée de Fourier discrète, calculée à l'aide d'un algorithme FFT, fournit directement les coefficients de Fourier (see here). Par contre, si vous avez une expression analytique pour la fonction, vous avez probablement besoin d'un solveur mathématique symbolique quelconque.

Answer 3

Le background information sur la routine est la première ressource à laquelle vous devriez vous référer.

Answer 4

Numpy n'est pas vraiment le bon outil pour calculer les composants de la série Fourier, car vos données doivent être échantillonnées discrètement. Vous voulez vraiment utiliser quelque chose comme Mathematica ou utiliser des transformées de Fourier. Pour le faire en gros, regardons quelque chose de simple une onde triangulaire de période 2pi, où l'on peut facilement calculer les coefficients de Fourier (c_n = -i ((-1)^(n + 1))/n pour n > 0, par exemple, c_n = {-i, i/2, -i/3, i/4, -i/5, i/6, ...} pour n = 1,2,3,4,5, 6 (en utilisant Somme (c_n exp (i 2 pi nx)) en série de Fourier)

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 

Si vous regardez les premières composantes de Fourier:.

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, 

négligence d'abord les composants qui sont proches 0 en raison de la précision en virgule flottante (~ 1e-16, étant zéro). La partie la plus difficile est de voir que les nombres 3.14159 (qui ont surgi avant que nous divisons par la période d'un 1000) devraient également être identifiés comme zéro, car la fonction est périodique). Donc, si nous négligeons ces deux facteurs que nous obtenons:

fourier_trans = [0,0,-i,0,i/2,0,-i/3,0,i/4,0,-i/5,0,-i/6, ... 

et vous pouvez voir les numéros de série de Fourier viennent comme tout autre nombre (je ne l'ai pas étudié, mais je crois que les composants correspondent à [c0, c- 1, c1, c-2, c2, ...]). J'utilise les conventions selon le wiki: http://en.wikipedia.org/wiki/Fourier_series.

Encore une fois, je suggère d'utiliser mathematica ou un système d'algèbre informatique capable d'intégrer et de traiter des fonctions continues.

Answer 5

Comme d'autres réponses l'ont mentionné, il semble que ce que vous cherchez est un package informatique symbolique, donc numpy ne convient pas. Si vous souhaitez utiliser une solution gratuite basée sur python, alors sympy ou sage devrait répondre à vos besoins.

Answer 6

Ceci est une vieille question, mais comme je devais coder ceci, je poste ici la solution qui utilise le module numpy.fft, qui est probablement plus rapide que d'autres solutions artisanales.

Le DFT est le bon outil pour le calcul de la précision numérique des coefficients de la série de Fourier d'une fonction, définie comme une expression analytique de l'argument ou comme une fonction d'interpolation numérique sur certains points discrets.

C'est la mise en œuvre, ce qui permet de calculer les coefficients à valeurs réelles de la série de Fourier, ou les coefficients valeurs complexes, en passant un approprié 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 

Ceci est un exemple d'utilisation:

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

Et l'intrigue des coefficients a0,a1,...,a10,b1,b2,...,b10 résultant:

Ceci est un test facultatif pour la fonction, pour les deux modes de fonctionnement. Vous devez exécuter ceci après l'exemple, ou définir une fonction périodique f et une période T avant d'exécuter le 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)