I've been working on understanding Fourier analysis both with pencil and paper as well as on how to write an algorithm that computes fourier coefficients (yes, I know there are already libraries for this sort of thing, that is not the point).
I have been reading about how symmetry in the function being represented can determine simplification of some of the coefficients. Even functions give us cosine Fourier series, odd ones give us sine series, and functions that are neither have contributions from both sine and cosine.
So I wrote my algorithm and graphed three function, even, odd and neither. At first everything looked great, but then I noticed something odd. When I extended the graph beyond the interval that I used to solve for my coefficients, that the symmetry matched that of the original function. However, it looks like every other period, the graph is reflected about the horizontal axis. Is this correct? I thought that the graph should repeat itself exactly. What is the reason behind this reflection?
I used the complex definition of Fourier series coefficients:
$$c_m = \frac{1}{2L}\sum_{-\infty}^{\infty}c_n\delta_{n,m}\int_{-L}^Ldx = \frac{1}{2L}\int_{-L}^L f(x) \exp(\frac{-im\pi x}{L})$$
To arrive at the following graphs.
$$x^2 \text{ even plot:} $$
$$x^3 : \text{odd plot}$$
$$e^x \text{ neither even nor odd: }$$
def coefficients(fn, dx, L):
"""
Calculate the complex form fourier series coefficients for the first M
waves using extended trapezoid rule.
: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 Fourier 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)
print(m)
coeffs = np.zeros(m, dtype=np.complex_)
xk = np.arange(-L, L, 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, 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):
n = iw - w/2
s[i] += c_coef[iw] * np.exp(1j * n * np.pi * ix / L)
# If a plot is desired:
if plot:
plt.suptitle("Fourier Series Plot")
plt.xlabel(r"$x$")
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__':
# step size for calculating c_m coefficients with trapezoid rule
deltax = .025 * np.pi
# length of interval for Fourier Series is 2*l
l = np.pi
c_m = coefficients(lambda x: np.e**x, deltax, l)
# The x range we would like to interpolate function values
x = np.arange(-l*4, l*4, .01)
sol = fourier_graph(x, l, c_m, lambda x: x**2, plot=True)
Best Answer
If you print the values generated by the line
n = iw - w/2, you will see that your script is generating non-integer values of $n$, specifically $n = -20.5, -19.5, \ldots, -0.5, 0.5, \ldots, 19.5$. Certainly nothing stops you from computing Fourier coefficients using these values: $$c_n = \frac{1}{2L}\int_{-L}^{L} f(x) \exp(-i \pi n x / L)\ dx$$ But when you attempt to recover the original function via $$f(x) = \sum_{n}c_n \exp(i \pi n x / L)$$ it doesn't work. The original function $f$ has period $2L$ (which in your case is $2\pi$), so it satisfies $f(x + 2L) = f(x)$. But plugging in $x+2L$ in place of $x$ on the right hand side in the above equation gives us $$f(x+2L) = \sum_{n}c_n \exp(i \pi n (x+2L) / L) = \sum_{n}c_n \exp(i \pi n x / L) \exp(i 2 \pi n )$$ If $n$ is an integer then that last factor is $\exp(i 2\pi n) = 1$ so you get the same result as for $f(x)$. But in your case the values of $n$ are equal to an integer plus $1/2$ (say $n = n' + 1/2$ where $n'$ is an integer) so you get $\exp(i 2\pi (n' + 1/2)) = \exp(i 2\pi n') \cdot \exp(i \pi) = 1 \cdot -1 = -1$, resulting in the incorrect result $f(x+2L) = -f(x)$.If you change your script so that you use integer values of $n$ (for example by changing
n = iw - w/2ton = iw - (w-1)/2) then this extra factor of $-1$ goes away and you get the expected result.Your script has some other problems, such as losing the imaginary part of the
sarray due to not usingdtype=np.complex_. But the main mathematical issue is the one explained above.