Implementation of Fourier Series up to ‘n’ Harmonics in MATLAB

In Mathematical Calculus, the expanded form of a periodic function “f(x)” in terms of an infinite sum of cosines & sines is called as Fourier series. It makes use of the orthogonality relationships of the cosine & sine functions. In other words, the Fourier series can also be defined as a way of representing any periodic function “f(x)” as a sum of cosine and sine functions (possibly infinite). 

For a periodic function “f(x)”, the Fourier series of “f(x)” in the Range of [ c, c+2l ] is given by :

 (In General, l=π)  

Here:

Here, a_n & b_n are called Fourier cosine and sine coefficients respectively.  

Note: If in the above formula of Fourier Series, instead of Infinity we use summation from n=1 to n=k then we call it as Fourier series of f(x) up to ‘k’ harmonics.

MATLAB functions used in the code are:

• disp(“txt”): This Method displays the Message-“txt” to the User.
• input(“txt”): This Method displays the Message-“txt” and waits for the user to input a value and press the Return key.
• int(y,x₁,x₂): This Method computes the definite integral of ‘y’ from x₁ to x₂.
• vpa(f,4): This Method evaluates each element of ‘f’ to at least 4-Significant digits.
• ezplot(y,[x₁,x₂]): This Method plots ‘y’ over the specified interval [x₁,x₂].
• title(‘GFG’): This method adds the specified title-“GFG” to a plot.
• legend(A,B,…): This Method creates a legend with descriptive labels for each plotted line.
• strcat(A,B): This Method horizontally concatenates the text in its input arguments. 
• char(f): This Method converts the input array ‘f’ to a character array.

Example:

Matlab

% MATLAB code for Implementation of
% Fourier Series up to 'n' Harmonics in MATLAB:
clear all
clc
disp("Implementation of Fourier Series
up to 'n' Harmonics in MATLAB | GeeksforGeeks")
syms x

f=input("Enter the function of x, whose fourier series is to be found:");
I=input("Enter the limits of Integration [a , b] :");
k=input("Enter the number of Harmonics:");

% Lower limit of Integration
a=I(1);

 % Upper limit of Integration
b=I(2);
l=(b-a)/2;
a0=(1/l)*(int(f,a,b));
Fx=a0/2;

% Calculating the nth Harmonic
for n=1:k
%To creates a new figure window
% using default property values
    figure;
    an(n)=(1/l)*(int(f*cos(n*pi*x/l),a,b));
    bn(n)=(1/l)*(int(f*sin(n*pi*x/l),a,b));
    Fx=Fx+((an(n))*cos(n*pi*x/l))+((bn(n))*sin(n*pi*x/l));

    % To evaluate Each element of
     % Fx to at least 4-Significant digits
    Fx=vpa(Fx,4);

    % To plot the curve Fx (Fourier series upto
    % nth Harmonic) in the given interval [a,b]
    ezplot(Fx,[a,b]);

    % To add a second line plot (Given Function/f) without deleting
    % the existing line plot (Fourier series upto nth Harmonic/Fx)
    hold on;
    ezplot(f,[a,b]);

    % To plot the curve f (Given Function)
    % in the given interval [a,b]
    title(["The Fourier series upto ", num2str(n)," Harmonics is:"]);

    % To create a legend with descriptive labels
    % for both of the plotted lines (f & Fx)
    legend("Fourier series","Given Function");

   % To set the hold state to off
    hold off;
end
disp(strcat('The Fourier series upto ', num2str(n),' Harmonics is:',char(Fx)))

Output:

Input: For a function  in the Range= [ -π, π ],
       Fourier Series of f(x) up to '3' Harmonics is:

Input: For a function f(x)=e-x in the Range= [ 0, 2π ],
       Fourier Series of f(x) up to '4' Harmonics is: