transformada serie fourier espectro python numpy fft

python - serie - sympy fourier



¿Cómo calcular una serie de Fourier en Numpy? (6)

¿Tiene una lista de muestras discretas de su función, o su función es discreta? Si es así, la Transformada Discreta de Fourier, calculada utilizando un algoritmo FFT, proporciona los coeficientes de Fourier directamente ( ver aquí ).

Por otro lado, si tiene una expresión analítica para la función, probablemente necesite algún tipo de resolución simbólica de matemáticas.

Tengo una función periódica del período T y me gustaría saber cómo obtener la lista de los coeficientes de Fourier. Intenté usar el módulo fft de numpy pero parece más dedicado a las transformadas de Fourier que a las series. Tal vez sea una falta de conocimiento matemático, pero no puedo ver cómo calcular los coeficientes de Fourier a partir de fft.

Ayuda y / o ejemplos apreciados.

Al final, lo más simple (calcular el coeficiente con una suma riemann) fue la forma más portátil / eficiente / robusta de resolver mi 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)

me da

Como han mencionado otras respuestas, parece que lo que está buscando es un paquete de cómputo simbólico, por lo que el número no es adecuado. Si desea utilizar una solución gratuita basada en Python, sympy o sage deben satisfacer sus necesidades.

Esta es una pregunta antigua, pero como tuve que codificar esto, estoy publicando aquí la solución que usa el módulo numpy.fft , que probablemente sea más rápido que otras soluciones hechas a mano.

La DFT es la herramienta adecuada para el trabajo de cálculo hasta la precisión numérica de los coeficientes de la serie de Fourier de una función, definida como una expresión analítica del argumento o como una función de interpolación numérica sobre algunos puntos discretos.

Esta es la implementación, que permite calcular los coeficientes de valores reales de la serie de Fourier, o los coeficientes de valores complejos, pasando un return_complex apropiado:

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

Este es un ejemplo de uso:

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

Y la gráfica de los coeficientes a0,a1,...,a10,b1,b2,...,b10 :

Esta es una prueba opcional para la función, para ambos modos de operación. Debe ejecutar esto después del ejemplo, o definir una función periódica f un período T antes de ejecutar el código.

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

Numpy no es la herramienta correcta para calcular los componentes de la serie de Fourier, ya que sus datos deben ser muestreados discretamente. Realmente desea utilizar algo como Mathematica o debería usar transformadas de Fourier.

Para hacerlo más o menos, veamos algo simple como una onda triangular del período 2pi, donde podemos calcular fácilmente los coeficientes de Fourier (c_n = -i ((-1) ^ (n + 1)) / n para n> 0; por ejemplo, , c_n = {-i, i / 2, -i / 3, i / 4, -i / 5, i / 6, ...} para n = 1,2,3,4,5,6 (usando Sum (c_n exp (i 2 pi nx)) como series 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 nos fijamos en los primeros componentes 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,

Primero descuide los componentes que están cerca de 0 debido a la precisión del punto flotante (~ 1e-16, como cero). La parte más difícil es ver que los números 3.14159 (que surgieron antes de dividirnos por el período de 1000) también deben reconocerse como cero, ya que la función es periódica). Así que si descuidamos esos dos factores obtenemos:

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

y puede ver que los números de las series de Fourier aparecen como cualquier otro número (no he investigado, pero creo que los componentes corresponden a [c0, c-1, c1, c-2, c2, ...]). Estoy usando convenciones de acuerdo con wiki: http://en.wikipedia.org/wiki/Fourier_series .

De nuevo, sugeriría utilizar mathica o un sistema de álgebra computacional capaz de integrar y tratar funciones continuas.