Understanding Fourier Transforms in abstract math terms

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:

With 200 terms:

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)

• 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. – robert bristow-johnson Nov 12 '19 at 4:35
• 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? – Marcus Müller Nov 12 '19 at 8:08
• @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. – rocksNwaves Nov 12 '19 at 15:32
• @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. – rocksNwaves Nov 12 '19 at 15:44
• @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/… – rocksNwaves Nov 12 '19 at 15:56

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.

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