Code covered by the BSD License

# Fourier series of square wave. Demo of Gibbs phenomenon with overshoot calculation

### अंकित अ. भुराने Ankit A. Bhurane (view profile)

23 Sep 2013 (Updated )

This file gives a simple demonstration of how a square wave can be approximated by Fourier series.

Fourier_Series_Ankit_A_Bhurane.m
```%% ========================================================================
%% Title:
%  1) Fourier series approximation of square wave.
%  2) Demonstration of Gibbs phenomenon (verification of Fig. 3.9 of [1])
%% Author:
%  Ankit A. Bhurane (ankit.bhurane@gmail.com)
%% Expression:
% The Fourier series approximation of a square wave signal existing between
% -Tau/2 to Tau/2 and period of T0 will have the form:
%
% Original signal to be approximated:
%                         :
%               __________:__________  A
%              |          :          |
%              |          :          |
%              |          :          |
%    __________|          :          |__________
% -T0/2     -Tau/2        0        Tau/2       T0/2
%                         :
%                         :
%
% Its Fourier series approximation:
%
%                  Inf
%                  ___
%         A*Tau   \     / sin(pi*n*Tau/T0)                    \
% r(t) = -------   |   | ------------------ exp^(j*n*2*pi/T0*t) |
%          T0     /___  \  (pi*n*Tau/T0)                      /
%               n = -Inf
%
% The left term inside summation are the Fourier series coeffs (Cn). The
% right term is the Fourier series kernel.
% Tau: range of square wave, T: period of the square wave,
% t: time variable, n: number of retained coefficients.
%
%% Observations:
% 1) As number of retained coefficients tends to infinity, the approximated
% signal value at the discontinuity converge to half the sum of values on
% either side.
% 2) Ripples does not decrease with increasing coefficients with
% approximately 9% overshoot.
% 3) Energy in the error between original and approximated signal, reduces
% as the number of retained coefficients are increased.
%
%% References:
% [1] Oppenheim, Willsky, Nawab, "Signals and Systems", PHI, Second edition
% [2] Dean K. Frederick and A. Bruce Carlson, "Fourier series" section in
%     Linear systems in communication and control
%% Copyright (c) 2013-2014 | Ankit A. Bhurane
%% =========================================================================

clc; clear all; close all;

% Specification
A = 1; % Peak-to-peak amplitude of square wave
Tau = 10; % Total range in which the square wave is defined (here -5 to 5)
T0 = 20; % Period (time of repeatation of square wave), here 10
C = 30; % Coefficients (sinusoids) to retain
N = 1001; % Number of points to consider
t = linspace(-(T0-Tau),(T0-Tau),N); % Time axis
X = zeros(1,N); X(t>=-Tau/2 & t<=Tau/2) = A; % Original signal
R = 0; % Initialize the approximated signal
k = -C:C; % Fourier coefficient number axis
f = zeros(1,2*C+1); % Fourier coefficient values

% Loop for plotting approximated signals for different retained coeffs.
for c = 0:C % Number of retained coefficients
for n = -c:c % Summation range (See equation above in comments)

% Sinc part of the Fourier coefficients calculated separately
if n~=0
Sinc = (sin(pi*n*Tau/T0)/((pi*n*Tau/T0))); % At n NOTEQUAL to 0
else
Sinc = 1; % At n EQUAL to 0
end
Cn = (A*Tau/T0)*Sinc; % Actual Fourier series coefficients
f(k==n) = Cn; % Put the Fourier coefficients at respective places
R = R + Cn*exp(1j*n*2*pi/T0.*t); % Sum all the coefficients
end

R = real(R); % So as to get rid of 0.000000000i (imaginary) factor
Max = max(R); Min = min(R); M = max(abs(Max),abs(Min)); % Maximum error
Overshoot = ((M-A)/A)*100; % Overshoot calculation
E = sum((X-R).^2); % Error energy calculation

% Plots:
% Plot the Fourier coefficients
subplot 121; stem(k,f,'m','LineWidth',3); axis tight; grid on;
xlabel('Fourier coefficient index');ylabel('Magnitude');
title('Fourier coefficients');

% Plot the approximated signal
subplot 122; plot(t,X,t,R,'m','LineWidth',3); axis tight; grid on;
xlabel('Time (t)'); ylabel('Amplitude');
title(['Approximation for N = ', num2str(c),...
'. Overshoot = ',num2str(Overshoot),'%','. Error energy: ',num2str(E)])

pause(0.1); % Pause for a while
R = 0; % Reset the approximation to calculate new one
end

%% ========================================================================```