Fourier Series of a given function

Question

This is a very simple code that expresses a function in terms of Trigonometric Fourier Series and generates a animation based on the harmonics. I would like to know some ways to improve the performance of the code and its readability. I believe that any question about the mathematical part can be answered here.

Here's the whole code:

import matplotlib.pyplot as plt
from numpy import sin, cos, pi, linspace
from scipy.integrate import quad as integral
from celluloid import Camera
def desired_function(t):
 return t**2
def a_zero(function, half_period, inferior_limit, superior_limit):
 return (1/(2*half_period)) * integral(function, inferior_limit, superior_limit)[0]
def a_k(k, half_period, inferior_limit, superior_limit):
 return (1/half_period) * integral(lambda x: desired_function(x) * cos(x * k * (pi / half_period)), inferior_limit, superior_limit)[0]
def b_k(k, half_period, inferior_limit, superior_limit):
 return (1/half_period) * integral(lambda x: desired_function(x) * sin(x * k * (pi / half_period)), inferior_limit, superior_limit)[0]
def main():
 fig = plt.figure()
 camera = Camera(fig)
 inferior_limit, superior_limit = -pi, pi
 half_period = (superior_limit - inferior_limit) / 2
 x = linspace(inferior_limit, superior_limit, 1000)
 f = a_zero(desired_function, half_period, inferior_limit, superior_limit)
 # The sum ranging from 1 to total_k
 total_k = 30
 for k in range(1, total_k+1):
 f += a_k(k, half_period, inferior_limit, superior_limit)*cos(k*x*(pi/half_period)) +\
 b_k(k, half_period, inferior_limit, superior_limit)*sin(k*x*(pi/half_period))
 plt.plot(x, desired_function(x), color='k')
 plt.plot(x, f, label=f'k = {k}')
 camera.snap()
 animation = camera.animate()
 plt.close()
 animation.save('animation.gif')
if __name__ == '__main__':
 main()

And my biggest issue is this for loop being used as a Σ, is there a more elegant way to express this?

# The sum ranging from 1 to total_k
 total_k = 30
 for k in range(1, total_k+1):
 f += a_k(k, half_period, inferior_limit, superior_limit)*cos(k*x*(pi/half_period)) +\
 b_k(k, half_period, inferior_limit, superior_limit)*sin(k*x*(pi/half_period))

Thanks in advance to anyone who takes the time to help me with this!

Answer 1

I like it. This code hews closely to the underlying math, using lovely identifiers, and is very clear.

---

For {a,b}_k, here's the only part I'm not super happy with:

def a_k(k, half_period, inferior_limit, superior_limit):
 return (1/half_period) * integral(lambda x: desired_function(x) * cos(x * k * (pi / half_period)), inferior_limit, superior_limit)[0]

We went with desired_function as a global. Ok, not the end of the world. Just sayin', consider adding it to the parameters.

I probably would have gone with a concise lo, hi instead of inferior_limit, superior_limit. But hey, that's just me, it's a matter of taste. Don't change it on my account, looks good as it is. Maybe package them both as the 2-tuple limits, and then the call goes from

 ... * integral(... * (pi / half_period)), inferior_limit, superior_limit)[0]

to

 ... * integral(... * (pi / half_period)), *limits)[0]

---

def main():
 ...
 x = linspace(inferior_limit, superior_limit, 1000)
 ...
 total_k = 30

We have two magic numbers, both quite reasonable.

Consider burying one or both of them in the signature:

def main(num_points=1000, total_k=30):

One benefit of such practice is it makes it easy for unit tests to adjust numbers to more convenient values.

---

It wouldn't hurt to sometimes add a """docstring""" to a function, explaining what it expects and what it offers. Not every function, mind you, just the big one(s).

---

To follow PEP-8 conventions more closely, consider running black over this source. Whitespace within math expressions tends to improve readability.

---

more elegant way to express [the summation?]

No, I respectfully disagree with your assessment. It exactly corresponds to the math and it looks beautiful to me. You could change it, but then I fear it would be less clear.

Now, if we wanted to worry about efficiency and execution speed, there are some constants that could be hoisted out of the loop, and the repeated calls to integral() seem like they might show up on cProfile's hit list. But if we cared about elapsed time we would use numba's @jit or another technique, and compilers are good at hoisting. This code achieves its design goals and is maintainable. ATM it looks good to me.

Ship it!

Answer 2

You have really clean code, I can only suggest a few improvements:

The functions a_k and b_k have a quite long line, I suggest using a named function rather than a lambda like this:

def a_k(k, half_period, inferior_limit, superior_limit):
 def integrand(x):
 return desired_function(x) * cos(x * k * (pi / half_period))
 return (1/half_period) * integral(integrand, inferior_limit, superior_limit)[0]

In this way you have shorter line of code and a self-documenting name following standard mathematical conventions for the integrand

The same for b_k.

---

Regarding the summmation in the loop I would rewrite it from:

total_k = 30
for k in range(1, total_k+1):
 f += a_k(k, half_period, inferior_limit, superior_limit)*cos(k*x*(pi/half_period)) +\
 b_k(k, half_period, inferior_limit, superior_limit)*sin(k*x*(pi/half_period))

To:

total_k = 30
for k in range(1, total_k+1):
 trig_argument = k*x*(pi/half_period)
 cosine_contribution = a_k(k, half_period, inferior_limit, superior_limit)*cos(trig_argument)
 sine_contribution = b_k(k, half_period, inferior_limit, superior_limit)*sin(trig_argument)
 f += cosine_contribution + sine_contribution

To add descriptive names to variables and highlight that sine and cosine are called with the same argument, but at the cost of being a bit longer. Anyhow also your version is fine. Another possibility is integration_arguments = (k, half_period, inferior_limit, superior_limit) and using it later, to highlight that the integration is over the same bounds for sine and cosine, this would also highlight the similarity of the operations for sine and cosine.