How to calculate a Fourier series in Numpy?

Question

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.

Answer 1

The background information about the routine is the first resource you should refer to.

Answer 2

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:

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:

Answer 3

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.

Answer 4

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:

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)

Answer 5

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.

Answer 6

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.