Я имею периодическую функцию периода T и хотел бы знать, как получить список коэффициентов Фурье. Я попытался использовать fft модуль из numpy, но, похоже, он больше посвящен преобразованиям Фурье, чем серии. Может быть, это недостаток математических знаний, но я не вижу, как вычислить коэффициенты Фурье из fft.
Справка и/или примеры, оцененные.
В конце концов, самая простая вещь (вычисление коэффициента с суммой римана) была самым переносимым/эффективным/надежным способом решения моей проблемы:
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)
дает мне:
Numpy - это не самый подходящий инструмент для расчета компонентов серии Fourier, так как ваши данные должны быть дискретно отобраны. Вы действительно хотите использовать что-то вроде Mathematica или должны использовать преобразования Фурье.
Чтобы грубо сделать это, рассмотрим что-то простую треугольную волну периода 2pi, где мы можем легко вычислить коэффициенты Фурье (c_n = -i ((-1) ^ (n + 1))/n при 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.
Опять же, я бы предложил использовать математику или систему компьютерной алгебры, способную интегрировать и работать с непрерывными функциями.
Как упоминалось в других ответах, кажется, что то, что вы ищете, является символическим вычислительным пакетом, поэтому numpy не подходит. Если вы хотите использовать бесплатное решение на основе python, то sympy или sage должен отвечать вашим потребностям.
У вас есть список дискретных образцов вашей функции, или ваша функция сама дискретна? Если это так, дискретное преобразование Фурье, рассчитанное с использованием алгоритма БПФ, непосредственно дает коэффициенты Фурье (см. Здесь).
С другой стороны, если у вас есть аналитическое выражение для функции, вам, вероятно, нужен какой-то символический математический решатель.
Это старый вопрос, но поскольку мне пришлось его кодировать, я размещаю здесь решение, которое использует модуль 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)
справочная информация о подпрограмме - это первый ресурс, на который вы должны ссылаться.