Come calcolare una serie di Fourier in Numpy?

Question

Ho una funzione periodica del periodo T e vorrei sapere come ottenere l'elenco dei coefficienti di Fourier. Ho provato a usare il modulo fft da numpy ma sembra più dedicato alle trasformate di Fourier che alle serie. Forse è una mancanza di conoscenza matematica, ma non riesco a vedere come calcolare i coefficienti di Fourier da fft.

Aiuto e/o esempi apprezzati.

Answer 1

Alla fine, la cosa più semplice (il calcolo del coefficiente con una somma di Riemann) è stato il/modo più portatile efficiente/robusto per risolvere il mio problema:

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) 

mi da:

Answer 2

Avete una lista di campioni discreti della vostra funzione o la vostra funzione è discreta? In tal caso, la trasformazione di Fourier discreta, calcolata utilizzando un algoritmo FFT, fornisce direttamente i coefficienti di Fourier (see here).

D'altra parte, se si dispone di un'espressione analitica per la funzione, è probabilmente necessario un simbolico risolutore matematico di qualche tipo.

Answer 3

Il background information sulla routine è la prima risorsa a cui è necessario fare riferimento.

Answer 4

Numpy non è lo strumento giusto per calcolare i componenti della serie di Fourier, in quanto i dati devono essere campionati in modo discreto. Vuoi davvero usare qualcosa come Mathematica o usare le trasformazioni di Fourier.

Per fare approssimativamente, diamo un'occhiata a qualcosa di semplice un'onda triangolare di punto 2pi, in cui possiamo facilmente calcolare i coefficienti di Fourier (c_n = -i ((-1)^(n + 1))/n per n > 0; ad esempio, c_n = {-i, i/2, -i/3, i/4, -i/5, i/6, ...} per n = 1,2,3,4,5, 6 (utilizzando Sum (c_n exp (i 2 pi NX)) come serie di 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 

Se si guardano le prime componenti di 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, 

Prima trascurare le componenti che sono vicino 0 a causa della precisione in virgola mobile (~ 1e-16, come zero) la parte più difficile è vedere che i numeri 3.14159 (che sono sorti prima di dividere per il periodo di 1000) dovrebbero essere riconosciuti come zero, poiché la funzione è periodica). Quindi, se trascuriamo questi due fattori si ottiene:

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

ed è possibile vedere i numeri serie di Fourier venire come ogni altro numero (non ho indagato, ma credo che i componenti corrispondano a [C0, c- 1, c1, c-2, c2, ...]). Sto usando le convenzioni secondo wiki: http://en.wikipedia.org/wiki/Fourier_series.

Ancora una volta, suggerirei di utilizzare la matematica o un sistema di algebra informatico in grado di integrare e gestire funzioni continue.

Answer 5

Come altre risposte hanno menzionato, sembra che quello che stai cercando sia un pacchetto di calcolo simbolico, quindi non è adatto a numpy. Se desideri utilizzare una soluzione basata su Python gratuita, allora sympy o sage dovrebbero soddisfare le tue esigenze.

Answer 6

Questa è una vecchia domanda, ma dal momento che ho dovuto codificare questo, sto postando qui la soluzione che utilizza il modulo numpy.fft, che è probabilmente più veloce di altre soluzioni artigianali.

Lo DFT è lo strumento giusto per calcolare fino alla precisione numerica i coefficienti della serie di Fourier di una funzione, definita come un'espressione analitica dell'argomento o come funzione di interpolazione numerica su alcuni punti discreti.

Questa è l'attuazione, che permette di calcolare i coefficienti a valori reali della serie di Fourier, oi valori complessi coefficienti, passando un appropriato 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 

Questo è un esempio di utilizzo:

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

E la trama delle conseguenti a0,a1,...,a10,b1,b2,...,b10 coefficienti:

Questo è un test opzionale per la funzione, per entrambe le modalità di funzionamento. È necessario eseguire questo dopo l'esempio oppure definire una funzione periodica f e un periodo T prima di eseguire il codice.

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