Fourier series of square wave. Demo of Gibbs phenomenon with overshoot calculation
23 Sep 2013
25 Sep 2013)
This file gives a simple demonstration of how a square wave can be approximated by Fourier series.
% 1) Fourier series approximation of square wave.
% 2) Demonstration of Gibbs phenomenon (verification of Fig. 3.9 of )
% Ankit A. Bhurane (firstname.lastname@example.org)
% 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:
% 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.
% 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.
%  Oppenheim, Willsky, Nawab, "Signals and Systems", PHI, Second edition
%  Dean K. Frederick and A. Bruce Carlson, "Fourier series" section in
% Linear systems in communication and control
%% Last Modified: Sept 24, 2013.
%% Copyright (c) 2013-2014 | Ankit A. Bhurane
clc; clear all; close all;
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
Sinc = (sin(pi*n*Tau/T0)/((pi*n*Tau/T0))); % At n NOTEQUAL to 0
Sinc = 1; % At n EQUAL to 0
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
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
% Plot the Fourier coefficients
subplot 121; stem(k,f,'m','LineWidth',3); axis tight; grid on;
xlabel('Fourier coefficient index');ylabel('Magnitude');
% 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