1
$\begingroup$

I am having a hard time implementing a method that computes Fourier transform coefficients for the complex form using the trapezoid rule.

I have floated questions in the math and stackoverflow networks, without getting any kind of answers that I can understand and implement.

Essentially, I am having the following issues:

  1. My DFT estimate of any given function seems to be at least off by half. I can correct this by multiplying each DFT coefficient by 2, but that's probably a terrible thing to do in lieu of finding the actual bug.
  2. No matter how many coefficients I calculate, the DFT does not converge very well to the exact function (unless I choose to estimate a periodic function).
  3. when I add more and more coefficients, things begin to look wrong after a certain point.

EDIT: Added two plots to illustrate issue #3.

With about 100 terms: enter image description here

With 200 terms: enter image description here

Unfortunately, you will have to keep the language of your help tailored towards someone who has only done abstract math. Here is the code:

import matplotlib.pyplot as plt
import numpy as np


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

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

    N = 2*L / dx

    # m is the number of DFT coefficients to calculate. N/2 coefficients are sufficient
    # Proof:
    # Minimum sampling given a oscillating function is lambda/2, where lambda is the
    # period of oscillation of the function we sample. In other words dx = lambda/2
    # for a minimally resolved  estimate. So with fixed dx, the smallest wavelength
    # we can hope to resolve is lambda_min = 2*dx. Also note that the wavelength of
    # sin(n*pi*x / L) is given by 2L/n. So then our minimum wavelength is given by
    # lambda_min = 2L/n_min = 2*dx. Rearranging, we see that the smallest wave number
    # n is n_min = 2L/lambda_min = 2L/(2*dx) = L/dx = L/(2L/N) = N/2.

    m = int(N/2 + 1)

    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] = 2/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 = 30

    # 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, 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)
$\endgroup$
  • 1
    $\begingroup$ i just don't wanna look at your code. and i don't know Python anyway. could you use $\LaTeX$ to express your mathematical problem in traditional mathematical notation? usually, factor of two off is due to the double-sided Fourier transform (with both positive and negative frequency) being related to a single-sided with real sinusoids having only positive frequencies. because: $$ \cos(\omega t) = \tfrac12 (e^{i \omega t} + e^{-i \omega t}) $$ sometimes a factor of 2 sneaks in there. $\endgroup$ – robert bristow-johnson Nov 12 at 4:35
  • $\begingroup$ Could you elaborate what you mean with "the DFT doesn't converge very well to the exact function": Which function? I presume you mean a function on $\mathbb R$, so what kind of measure are you applying to define convergence "quality" (I'm not familiar with "well-converging" as a term, but my analysis education lies a while in the past, so maybe that's a well-defined math term) for a finite sequence of complex numbers and a function on the reals? What is the "certain point" after which things look wrong? $\endgroup$ – Marcus Müller Nov 12 at 8:08
  • $\begingroup$ @MarcusMüller Sure! Well basically the relative error remains high even as more and more terms are added to the finite sequence. I was under the impression that for any smooth function, the Fourier Series should converge to to the exact function as the number of terms approaches infinity. This might have been a misconception on my part. Someone mentioned to me yesterday that functions that are not "square-integrable" will not necessarily converge. I haven't looked up what that means yet. $\endgroup$ – rocksNwaves Nov 12 at 15:32
  • $\begingroup$ @MarcusMüller and to answer your second question, given that my step size for calculating the coefficients is $0.025 \pi$, my plot starts to look wack when I exceed about 100 terms. By the time I get to 200, terms it's totally crazy. I'll post the plots in the question. $\endgroup$ – rocksNwaves Nov 12 at 15:44
  • $\begingroup$ @robertbristow-johnson Sure thing, here is a link to the question I posted on the math network, although I believe at this point the question about math has been answered and now my question is one of implementation: math.stackexchange.com/questions/3420961/… $\endgroup$ – rocksNwaves Nov 12 at 15:56
1
$\begingroup$

Big thanks to @user753642 for spotting my mistakes over on the stackoverflow network:

I was computing the $c_n$ coefficients from $n=0 \dots m$, where m is the number of wave functions in the sum. But by definitions the coefficients look like:

$$c_m = \frac{1}{2L}\sum_{n = -\infty}^{\infty}c_n\delta_{n,m}\int_{-L}^Ldx = \frac{1}{2L}\int_{-L}^L f(x) \exp(\frac{-im\pi x}{L})$$

And the reconstructed function looks like:

$$f(x) = \sum_{n=-\infty}^{\infty}c_n \exp(\frac{in \pi x}{L}) $$

So when I was reconstructing my function with the Fourier Series with $c_n$ with $n=0…m$, I should have been going from $n=−m/2…m/2$.

This fixed everything, and allowed me to drop the weird $N/2$ factor that I was using. This was achieved by changing the above code int two places. First when calculating the coefficients:

# 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

And then again when reconstructing the function:

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 * np.pi * ix / L)

I have a lot to learn about Fourier Series/Transforms. Thanks for all the helpful comments, everyone.

$\endgroup$
  • 1
    $\begingroup$ You can accept your own answer if you like. $\endgroup$ – Matt L. Nov 12 at 18:50
  • $\begingroup$ Apparently I have to weight 24 hours. $\endgroup$ – rocksNwaves Nov 13 at 15:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.