Calculate the Fourier series with a trigonometric approach

Question

I am trying to implement the function of Fourier series according to the following formulas:

... where ...

... and ...

Here is my approach to the problem:

import numpy as np import pylab as py # Define "x" range. x = np.linspace(0, 10, 1000) # Define "T", ie functions' period. T = 2 L = T / 2 # "f(x)" function definition. def f(x): return np.sin(np.pi * 1000 * x) # "a" coefficient calculation. def a(n, L, accuracy = 1000): a, b = -L, L dx = (b - a) / accuracy integration = 0 for i in np.linspace(a, b, accuracy): x = a + i * dx integration += f(x) * np.cos((n * np.pi * x) / L) integration *= dx return (1 / L) * integration # "b" coefficient calculation. def b(n, L, accuracy = 1000): a, b = -L, L dx = (b - a) / accuracy integration = 0 for i in np.linspace(a, b, accuracy): x = a + i * dx integration += f(x) * np.sin((n * np.pi * x) / L) integration *= dx return (1 / L) * integration # Fourier series. def Sf(x, L, n = 10): a0 = a(0, L) sum = 0 for i in np.arange(1, n + 1): sum += ((a(i, L) * np.cos(n * np.pi * x)) + (b(i, L) * np.sin(n * np.pi * x))) return (a0 / 2) + sum # x axis. py.plot(x, np.zeros(np.size(x)), color = 'black') # y axis. py.plot(np.zeros(np.size(x)), x, color = 'black') # Original signal. py.plot(x, f(x), linewidth = 1.5, label = 'Signal') # Approximation signal (Fourier series coefficients). py.plot(x, Sf(x, L), color = 'red', linewidth = 1.5, label = 'Fourier series') # Specify x and y axes limits. py.xlim([0, 10]) py.ylim([-2, 2]) py.legend(loc = 'upper right', fontsize = '10') py.show() 

... and this is what I get after building the result:

I read How to calculate the Fourier series in Numpy? , and I have already implemented this approach. It works fine, but it uses the exponentiality method, where I want to focus on trigonometry functions and the rectangular method when calculating the integrals for the coefficients a_{n} and b_{n} .

Thanks in advance.

UPDATE (SOLVED)

Finally, here is a working code example. However, I will spend more time on it, so if there is something that can be improved, it will be done.

 from __future__ import division import numpy as np import pylab as py # Define "x" range. x = np.linspace(0, 10, 1000) # Define "T", ie functions' period. T = 2 L = T / 2 # "f(x)" function definition. def f(x): return np.sin((np.pi) * x) + np.sin((2 * np.pi) * x) + np.sin((5 * np.pi) * x) # "a" coefficient calculation. def a(n, L, accuracy = 1000): a, b = -L, L dx = (b - a) / accuracy integration = 0 for x in np.linspace(a, b, accuracy): integration += f(x) * np.cos((n * np.pi * x) / L) integration *= dx return (1 / L) * integration # "b" coefficient calculation. def b(n, L, accuracy = 1000): a, b = -L, L dx = (b - a) / accuracy integration = 0 for x in np.linspace(a, b, accuracy): integration += f(x) * np.sin((n * np.pi * x) / L) integration *= dx return (1 / L) * integration # Fourier series. def Sf(x, L, n = 10): a0 = a(0, L) sum = np.zeros(np.size(x)) for i in np.arange(1, n + 1): sum += ((a(i, L) * np.cos((i * np.pi * x) / L)) + (b(i, L) * np.sin((i * np.pi * x) / L))) return (a0 / 2) + sum # x axis. py.plot(x, np.zeros(np.size(x)), color = 'black') # y axis. py.plot(np.zeros(np.size(x)), x, color = 'black') # Original signal. py.plot(x, f(x), linewidth = 1.5, label = 'Signal') # Approximation signal (Fourier series coefficients). py.plot(x, Sf(x, L), '.', color = 'red', linewidth = 1.5, label = 'Fourier series') # Specify x and y axes limits. py.xlim([0, 5]) py.ylim([-2.2, 2.2]) py.legend(loc = 'upper right', fontsize = '10') py.show() 

Answer 1

Consider developing your code differently, in blocks. You should be surprised if such code worked on the first try. Debugging is one option, as @ tom10 said. Another option is to quickly prototype the code step by step in the interpreter, even better with ipython.

You expect b_1000 to be nonzero, since the input f(x) is a sine wave with 1000 in it. Do you also expect all other coefficients to be zero?

Then you should focus only on the function b(n, L, accuracy = 1000) . Looking at it, 3 things go wrong. Here are some suggestions.

• multiplication dx is inside the loop. Sure what?
• in the loop, should i be an integer on the right? Is this an integer? by prototyping or debugging you will discover this.
• be careful when you write (1/L) or a similar expression. If you use python2.7 you are mistaken. If not, at least use from __future__ import division at the top of your source. Read this PEP if you don't know what I'm talking about.

If you solve these three points, b() will work. Then think of a same way.