Fourier series Coefficients and Visualization [ Python Program ]

Today I wrote a code that calculates the Fourier Coefficients.

In case you don’t know what a Fourier Series is, then, basically it is a way of approximating or representing a periodic function by a series of simple harmonic(sine and cosine) functions.

You can check out the wikipedia for  more information on it: https://en.wikipedia.org/wiki/Fourier_series

For a waveform
f(x) with period 2L=2π/k, we have that k=2π/2L=π/L and nkx =nπx/L , Then we have
fourier series coefficients

So, here’s how my code works:
I created a function ‘Fourier’ that takes in 4 arguments.
First two arguments specify the range for which the original function that is to be approximated is defined for. These are indicated as li (initial limit),lf (final limit).

The third argument ‘n’ is the number of harmonic terms to be used. (Larger the ‘n’ => better the approximation.)
And the final argument is the function ‘f’ which is to be approximated.

The function returns the Fourier coefficients based on formula shown in the above image. The coefficients are returned as a python list: [a0/2,An,Bn].
a0/2 is the first Fourier coefficient and is a scalar.
An and Bn are numpy 1d arrays of size n, which store the coefficients of cosine and sine terms respectively.

To calculate these coefficients I perform integration using the script.integrate module.

CODE

#Fourier Series Coefficients
#The following function returns the fourier coefficients,'a0/2', 'An' & 'Bn'
#
#User needs to provide the following arguments:
#
#l=periodicity of the function f which is to be approximated by Fourier Series
#n=no. of Fourier Coefficients you want to calculate
#f=function which is to be approximated by Fourier Series
#
#*Some necessary guidelines for defining f:
#*The program integrates the function f from -l to l so make sure you define the function f correctly in the interval -l to l.
#
#for more information on Fourier Series visit: https://en.wikipedia.org/wiki/Fourier_series
#
#Written by: Manas Sharma([email protected])
#For more useful toolboxes and tutorials on Python visit: https://www.bragitoff.com/category/compu-geek/python/
def fourier(li, lf, n, f):
    l = (lf-li)/2 
    # Constant term
    a0=1/l*integrate.quad(lambda x: f(x), li, lf)[0]
    # Cosine coefficents
    A = np.zeros((n))
    # Sine coefficents
    B = np.zeros((n))
    
    for i in range(1,n+1):
        A[i-1]=1/l*integrate.quad(lambda x: f(x)*np.cos(i*np.pi*x/l), li, lf)[0]
        B[i-1]=1/l* integrate.quad(lambda x: f(x)*np.sin(i*np.pi*x/l), li, lf)[0]

    return [a0/2.0, A, B]

Now, let’s have a look at some real world examples. The following program can be used to approximate 4 popular periodic functions using the Fourier series. Just copy and paste it to run and play with it yourself.

# Author: Manas Sharma
# Website: www.bragitoff.com
# Email: [email protected]
# License: MIT

import numpy as np
import scipy.integrate as integrate






# Non-periodic sawtooth function defined for a range [-l,l]
def sawtooth(x):
    return x


# Non-periodic square wave function defined for a range [-l,l]
def square(x):
    if x>0:
        return np.pi
    else:
        return -np.pi



# Non-periodic triangle wave function defined for a range [-l,l]
def triangle(x):
    if x>0:
        return x
    else:
        return -x



# Non-periodic cycloid wave function defined for a range [-l,l]
def cycloid(x):
    return np.sqrt(np.pi**2-x**2)

#Fourier Series Coefficients
#The following function returns the fourier coefficients,'a0/2', 'An' & 'Bn'
#
#User needs to provide the following arguments:
#
#l=periodicity of the function f which is to be approximated by Fourier Series
#n=no. of Fourier Coefficients you want to calculate
#f=function which is to be approximated by Fourier Series
#
#*Some necessary guidelines for defining f:
#*The program integrates the function f from -l to l so make sure you define the function f correctly in the interval -l to l.
#
#for more information on Fourier Series visit: https://en.wikipedia.org/wiki/Fourier_series
#
#Written by: Manas Sharma([email protected])
#For more useful tutorials on Python visit: https://www.bragitoff.com/category/compu-geek/python/
def fourier(li, lf, n, f):
    l = (lf-li)/2 
    # Constant term
    a0=1/l*integrate.quad(lambda x: f(x), li, lf)[0]
    # Cosine coefficents
    A = np.zeros((n))
    # Sine coefficents
    B = np.zeros((n))
    
    for i in range(1,n+1):
        A[i-1]=1/l*integrate.quad(lambda x: f(x)*np.cos(i*np.pi*x/l), li, lf)[0]
        B[i-1]=1/l* integrate.quad(lambda x: f(x)*np.sin(i*np.pi*x/l), li, lf)[0]

    return [a0/2.0, A, B]

if __name__ == "__main__":


    # Limits for the functions
    li = -np.pi
    lf = np.pi

    # Number of harmonic terms
    n = 3

    # Fourier coeffficients for various functions
    coeffs = fourier(li,lf,n,sawtooth)
    print('Fourier coefficients for the Sawtooth wave\n')
    print('a0 ='+str(coeffs[0]))
    print('an ='+str(coeffs[1]))
    print('bn ='+str(coeffs[2]))
    print('-----------------------\n\n')

    coeffs = fourier(li,lf,n,square)
    print('Fourier coefficients for the Square wave\n')
    print('a0 ='+str(coeffs[0]))
    print('an ='+str(coeffs[1]))
    print('bn ='+str(coeffs[2]))
    print('-----------------------\n\n')

    coeffs = fourier(li,lf,n,triangle)
    print('Fourier coefficients for the Triangular wave\n')
    print('a0 ='+str(coeffs[0]))
    print('an ='+str(coeffs[1]))
    print('bn ='+str(coeffs[2]))
    print('-----------------------\n\n')

    coeffs = fourier(li,lf,n,cycloid)
    print('Fourier coefficients for the Cycloid wave\n')
    print('a0 ='+str(coeffs[0]))
    print('an ='+str(coeffs[1]))
    print('bn ='+str(coeffs[2]))
    print('-----------------------\n\n')

OUTPUT:

>>   python3 fourier.py
Fourier coefficients for the Sawtooth wave

a0 =0.0
an =[0. 0. 0.]
bn =[ 2.         -1.          0.66666667]
-----------------------


Fourier coefficients for the Square wave

a0 =0.0
an =[0. 0. 0.]
bn =[4.00000000e+00 3.11674695e-16 1.33333333e+00]
-----------------------


Fourier coefficients for the Triangular wave

a0 =1.5707963267948966
an =[-1.27323954e+00  4.99600361e-16 -1.41471061e-01]
bn =[0. 0. 0.]
-----------------------


Fourier coefficients for the Cycloid wave

a0 =2.4674011002723417
an =[ 0.89414547 -0.3336097   0.1850662 ]
bn =[0. 0. 0.]
-----------------------

It is quite evident, how the coefficients of sine terms are 0 if the function is even (triangular wave, cycloid), and the coefficients of cosine terms are 0 if the function is odd (square wave, sawtooth wave).

Now, let us try to visualise these results with the help of an animation created using the following Python code:

ANIMATION CODE:

# Author: Manas Sharma
# Website: www.bragitoff.com
# Email: [email protected]
# License: MIT

import numpy as np
import matplotlib.pyplot as plt 
from matplotlib.pyplot import *
import scipy.integrate as integrate

fig = figure(figsize=(7, 7), dpi=120)



# Function that will convert any given function 'f' defined in a given range '[li,lf]' to a periodic function of period 'lf-li' 
def periodicf(li,lf,f,x):
    if x>=li and x<=lf :
        return f(x)
    elif x>lf:
        x_new=x-(lf-li)
        return periodicf(li,lf,f,x_new)
    elif x<(li):
        x_new=x+(lf-li)
        return periodicf(li,lf,f,x_new)

# The periodic version of sawtooth function 
def sawtoothP(li,lf,x):
    return periodicf(li,lf,sawtooth,x)

# Non-periodic sawtooth function defined for a range [-l,l]
def sawtooth(x):
    return x

# The periodic version of square function
def squareP(li,lf,x):
    return periodicf(li,lf,square,x)

# Non-periodic square wave function defined for a range [-l,l]
def square(x):
    if x>0:
        return np.pi
    else:
        return -np.pi

# The periodic version of triangle function
def triangleP(li,lf,x):
    return periodicf(li,lf,triangle,x)

# Non-periodic triangle wave function defined for a range [-l,l]
def triangle(x):
    if x>0:
        return x
    else:
        return -x

# The periodic version of cycloid function
def cycloidP(li,lf,x):
    return periodicf(li,lf,cycloid,x)

# Non-periodic cycloid wave function defined for a range [-l,l]
def cycloid(x):
    return np.sqrt(np.pi**2-x**2)

#Fourier Series Coefficients
#The following function returns the fourier coefficients,'a0/2', 'An' & 'Bn'
#
#User needs to provide the following arguments:
#
#li,lf = Range over which the original function f which is to be approximated by Fourier Series is defined. Period is assumed to be lf-li
#n=no. of Fourier Coefficients you want to calculate
#f=function which is to be approximated by Fourier Series
#
#*Some necessary guidelines for defining f:
#*The program integrates the function f from li to lf so make sure you define the function f correctly in the interval li to lf.
#
#for more information on Fourier Series visit: https://en.wikipedia.org/wiki/Fourier_series
#
#Written by: Manas Sharma([email protected])
#For more useful toolboxes and tutorials on Python visit: https://www.bragitoff.com/category/compu-geek/python/
def fourierCoeffs(li, lf, n, f):
    l = (lf-li)/2 
    # Constant term
    a0=1/l*integrate.quad(lambda x: f(x), li, lf)[0]
    # Cosine coefficents
    A = np.zeros((n))
    # Sine coefficents
    B = np.zeros((n))
    
    for i in range(1,n+1):
        A[i-1]=1/l*integrate.quad(lambda x: f(x)*np.cos(i*np.pi*x/l), li, lf)[0]
        B[i-1]=1/l*integrate.quad(lambda x: f(x)*np.sin(i*np.pi*x/l), li, lf)[0]

    return [a0/2.0, A, B]

# This functions returns the value of the Fourier series for a given value of x given the already calculated Fourier coefficients
def fourierSeries(coeffs,x,l,n):
    value = coeffs[0]
    for i in range(1,n+1):
        value = value + coeffs[1][i-1]*np.cos(i*np.pi*x/l) +  coeffs[2][i-1]*np.sin(i*np.pi*x/l)
    return value
    



if __name__ == "__main__":

    # plt.style.use('dark_background')
    plt.style.use('seaborn')
    

    # Limits for the functions
    li = -np.pi
    lf = np.pi
    l = (lf-li)/2.0

    # Number of harmonic terms
    n = 1
    for n in range(1,10):

        plt.title('Fourier Series Approximation\nSawtooth Wave\n n = '+str(n))
        # plt.title('Fourier Series Approximation\nSquare Wave\n n = '+str(n))
        # plt.title('Fourier Series Approximation\nTriangular Wave\n n = '+str(n))
        # plt.title('Fourier Series Approximation\nCycloid\n n = '+str(n))
        
        # Fourier coeffficients for various functions
        coeffsSawtooth = fourierCoeffs(li,lf,n,sawtooth)
        coeffsTriangle = fourierCoeffs(li,lf,n,triangle)
        coeffsSquare = fourierCoeffs(li,lf,n,square)
        coeffsCycloid = fourierCoeffs(li,lf,n,cycloid)

    

    


        # Step size for plotting
        step_size = 0.05

        # Limits for plotting
        x_l = -np.pi*2
        x_u = np.pi*2

        # Sample values of x for plotting
        x = np.arange(x_l,x_u,step_size)
        y1 = [sawtoothP(li,lf,xi) for xi in x]
        y1_fourier = [fourierSeries(coeffsSawtooth,xi,l,n) for xi in x]
        y2 = [squareP(li,lf,xi) for xi in x]
        y2_fourier = [fourierSeries(coeffsSquare,xi,l,n) for xi in x]
        y3 = [triangleP(li,lf,xi) for xi in x]
        y3_fourier = [fourierSeries(coeffsTriangle,xi,l,n) for xi in x]
        y4 = [cycloidP(li,lf,xi) for xi in x]
        y4_fourier = [fourierSeries(coeffsCycloid,xi,l,n) for xi in x]

        x_plot =[]
        # Sawtooth
        y_plot1 = []
        y_plot1_fourier = []
        # Square
        y_plot2 = []
        y_plot2_fourier = []
        # Triangle
        y_plot3 = []
        y_plot3_fourier = []
        # Cycloid
        y_plot4 = []
        y_plot4_fourier = []
    
        x_l_plot = x_l - 13
        x_u_plot = x_l_plot + 14
        plt.xlim(x_l_plot,x_u_plot)
        plt.ylim(-6,7)

        for i in range(x.size):
            
            x_plot.append(x[i])
            # Actual function values
            y_plot1.append(y1[i])
            y_plot2.append(y2[i])
            y_plot3.append(y3[i])
            y_plot4.append(y4[i])
            # Values from fourier series
            y_plot1_fourier.append(y1_fourier[i])
            y_plot2_fourier.append(y2_fourier[i])
            y_plot3_fourier.append(y3_fourier[i])
            y_plot4_fourier.append(y4_fourier[i])

            #Sawtooth
            plt.plot(x_plot,y_plot1,c='darkkhaki',label='Sawtooth Wave')
            plt.plot(x_plot,y_plot1_fourier,c='forestgreen',label='Fourier Approximation')
            #Square
            #plt.plot(x_plot,y_plot2,c='tomato',label='Square Wave')
            #plt.plot(x_plot,y_plot2_fourier,c='maroon',label='Fourier Approximation')
            #Triangular
            #plt.plot(x_plot,y_plot3,c='orange',label = 'Triangular Wave')
            #plt.plot(x_plot,y_plot3_fourier,c='darkgoldenrod',label='Fourier Approximation')
            #Cycloid
            # plt.plot(x_plot,y_plot4,c='slateblue',label='Cycloid')
            # plt.plot(x_plot,y_plot4_fourier,c='teal',label='Fourier Approximation')
            
            x_l_plot = x_l_plot + step_size
            x_u_plot = x_u_plot + step_size
            plt.xlim(x_l_plot,x_u_plot)
            plt.pause(0.001)
            if i==0:
                plt.legend()

        
        plt.clf()
    
    plt.show()

NOTE: You will have to comment and de-comment some parts to visualise the different functions.

VISUALIZATION:

That’s it. I hope it was not too difficult to understand. If you have any questions I will be glad to answer them.

Periodic_identity_function

[wpedon id="7041" align="center"]

7 thoughts on “Fourier series Coefficients and Visualization [ Python Program ]

    1. Thanks a lot Ma’am!!!! I’m grateful for your encouragement!

  1. Respected sir
    It’s a masterpiece. Thank you for sharing this excellent works
    With Regards,
    Dr. Anathnath Ghosh

  2. Hi Sir.

    I am trying to visualize the fourier series that approximate a top hat function but the approximation function does a zig-zag between certain values which do not approximate the function.

  3. It was great to see a visualization of Fourier Transform.
    I did the same job, for image compression, in MATLAB software. I think doing Visualization for stock price or image compression makes it better for pedagogical uses.

  4. great and excellent !
    how this could be adapted to work on a time serie data instead of the function f ?

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.