Finding Fourier coefficients algorithm

Question

Ok, so I have been trying to code a "naive" method to calculate the coefficients for a standard Fourier Series in complex form. I am getting very close, I think, but there are some odd behaviors. This may be more of a math question than programming one, but I already asked on math.stackexchange and got zero answers. Here is my working code:

import matplotlib.pyplot as plt
import numpy as np


def coefficients(fn, dx, m, L):
    """
    Calculate the complex form fourier series coefficients for the first M
    waves.

    :param fn: function to sample
    :param dx: sampling frequency
    :param m: number of waves to compute
    :param L: We are solving on the interval [-L, L]
    :return: an array containing M Fourier coefficients c_m
    """

    N = 2*L / dx
    coeffs = np.zeros(m, dtype=np.complex_)
    xk = np.arange(-L, L + dx, dx)

    # Calculate the coefficients for each wave
    for mi in range(m):
        coeffs[mi] = 1/N * sum(fn(xk)*np.exp(-1j * mi * np.pi * xk / L))

    return coeffs


def fourier_graph(range, L, c_coef, function=None, plot=True, err_plot=False):
    """
    Given a range to plot and an array of complex fourier series coefficients,
    this function plots the representation.


    :param range: the x-axis values to plot
    :param c_coef: the complex fourier coefficients, calculated by coefficients()
    :param plot: Default True. Plot the fourier representation
    :param function: For calculating relative error, provide function definition
    :param err_plot: relative error plotted. requires a function to compare solution to
    :return: the fourier series values for the given range
    """
    # Number of coefficients to sum over
    w = len(c_coef)

    # Initialize solution array
    s = np.zeros(len(range))
    for i, ix in enumerate(range):
        for iw in np.arange(w):
            s[i] += c_coef[iw] * np.exp(1j * iw * np.pi * ix / L)

    # If a plot is desired:
    if plot:
        plt.suptitle("Fourier Series Plot")
        plt.xlabel(r"$t$")
        plt.ylabel(r"$f(x)$")
        plt.plot(range, s, label="Fourier Series")

        if err_plot:
            plt.plot(range, function(range), label="Actual Solution")
            plt.legend()

        plt.show()

    # If error plot is desired:
    if err_plot:
        err = abs(function(range) - s) / function(range)
        plt.suptitle("Plot of Relative Error")
        plt.xlabel("Steps")
        plt.ylabel("Relative Error")
        plt.plot(range, err)
        plt.show()

    return s


if __name__ == '__main__':

    # Assuming the interval [-l, l] apply discrete fourier transform:

    # number of waves to sum
    wvs = 50

    # step size for calculating c_m coefficients (trap rule)
    deltax = .025 * np.pi

    # length of interval for Fourier Series is 2*l
    l = 2 * np.pi

    c_m = coefficients(np.exp, deltax, wvs, l)

    # The x range we would like to interpolate function values
    x = np.arange(-l, l, .01)
    sol = fourier_graph(x, l, c_m, np.exp, err_plot=True)

Now, there is a factor of 2/N multiplying each coefficient. However, I have a derivation of this sum in my professor's typed notes that does not include this factor of 2/N. When I derived the form myself, I arrived at a formula with a factor of 1/N that did not cancel no matter what tricks I tried. I asked over at math.stackexchange what was going on, but got no answers.

What I did notice is that adding the 1/N decreased the difference between the actual solution and the fourier series by a massive amount, but it's still not right. so I tried 2/N and got even better results. I am really trying to figure this out so I can write a nice, clean algorithm for basic fourier series before I try to learn about Fast Fourier Transforms.

So what am I doing wrong here?

Answer 1

assuming c_n is given by A_n as in mathworld

idem c_n = 1/T \int_{-T/2}^{T/2}f(x)e^{-2ipinx/T}dx

we can compute (trivially) the coeffs c_n analytically (which is a good way to compare to your trapezoidal integral)

k = (1-2in)/2
c_n = 1/(4*pi*k)*(e^{2pik} - e^{-2pik})

So your coeffs are likely to be properly computed (the both wrong curves look alike)

Now notice that when you reconstitue f, you add the coeff c_0 up to c_m

But the reconstruction should occur with c_{-m} to c_m

So you are missing half of the coeffs.

Below a fix with your adaptated coefficients function and the theoretical coeffs

import matplotlib.pyplot as plt
import numpy as np


def coefficients(fn, dx, m, L):
    """
    Calculate the complex form fourier series coefficients for the first M
    waves.

    :param fn: function to sample
    :param dx: sampling frequency
    :param m: number of waves to compute
    :param L: We are solving on the interval [-L, L]
    :return: an array containing M Fourier coefficients c_m
    """

    N = 2*L / dx
    coeffs = np.zeros(m, dtype=np.complex_)
    xk = np.arange(-L, L + dx, dx)

    # Calculate the coefficients for each wave
    for mi in range(m):
        n = mi - m/2
        coeffs[mi] = 1/N * sum(fn(xk)*np.exp(-1j * n * np.pi * xk / L))

    return coeffs


def fourier_graph(range, L, c_coef, ref, function=None, plot=True, err_plot=False):
    """
    Given a range to plot and an array of complex fourier series coefficients,
    this function plots the representation.


    :param range: the x-axis values to plot
    :param c_coef: the complex fourier coefficients, calculated by coefficients()
    :param plot: Default True. Plot the fourier representation
    :param function: For calculating relative error, provide function definition
    :param err_plot: relative error plotted. requires a function to compare solution to
    :return: the fourier series values for the given range
    """
    # Number of coefficients to sum over
    w = len(c_coef)

    # Initialize solution array
    s = np.zeros(len(range), dtype=complex)
    t = np.zeros(len(range), dtype=complex)

    for i, ix in enumerate(range):
        for iw in np.arange(w):
            n = iw - w/2
            s[i] += c_coef[iw] * (np.exp(1j * n * ix * 2 * np.pi / L))
            t[i] += ref[iw] * (np.exp(1j * n * ix * 2 * np.pi / L))

    # If a plot is desired:
    if plot:
        plt.suptitle("Fourier Series Plot")
        plt.xlabel(r"$t$")
        plt.ylabel(r"$f(x)$")
        plt.plot(range, s, label="Fourier Series")

        plt.plot(range, t, label="expected Solution")
        plt.legend()

        if err_plot:
            plt.plot(range, function(range), label="Actual Solution")
            plt.legend()

        plt.show()

    return s

def ref_coefficients(m):
    """
    Calculate the complex form fourier series coefficients for the first M
    waves.

    :param fn: function to sample
    :param dx: sampling frequency
    :param m: number of waves to compute
    :param L: We are solving on the interval [-L, L]
    :return: an array containing M Fourier coefficients c_m
    """

    coeffs = np.zeros(m, dtype=np.complex_)

    # Calculate the coefficients for each wave
    for iw in range(m):
        n = iw - m/2
        k = (1-(1j*n)/2)
        coeffs[iw] = 1/(4*np.pi*k)* (np.exp(2*np.pi*k) - np.exp(-2*np.pi*k))
    return coeffs

if __name__ == '__main__':

    # Assuming the interval [-l, l] apply discrete fourier transform:

    # number of waves to sum
    wvs = 50

    # step size for calculating c_m coefficients (trap rule)
    deltax = .025 * np.pi

    # length of interval for Fourier Series is 2*l
    l = 2 * np.pi

    c_m = coefficients(np.exp, deltax, wvs, l)

    # The x range we would like to interpolate function values
    x = np.arange(-l, l, .01)
    ref = ref_coefficients(wvs)
    sol = fourier_graph(x, 2*l, c_m, ref, np.exp, err_plot=True)