square_wave_fourier_series_demo

%===========================================================================
% Square Wave Fourier Series Demo  (square_wave_fourier_series_demo.m)
%==========================================================================
% The user can design various square wave by determining its period, pulse 
% width, time shift, dc value, etc. Then the program can automatically
% compute its Fourier series representation, and plot its amplitude spectrum 
% and phase spectrum.
% 
% 
% 
% By Jing Tian Email: [e-mail]
% 
% 

function square_wave_fourier_series_demo

clc; close all; clear all;

%parameter of input square wave
T0 = 10;    % period
tau = 5;    % pulse width   
A = 1;
dc = 1;     % dc level
ts = 0;     % time shift, positive: move right. negative: move left
M = 4;      % How many period are shown

Nf = M*T0/tau;  %number of FS components, i.e., C(-Nf),...,C(-1),C(0),C(1),...C(Nf).
d = tau/T0;     % duty cycle
w0 = 2*pi/T0;   % frequency

%Figure 1, plot square wave
%The square signal is 'on' in the intervals [-tau/2,tau/2],
%while 'off' in the rest intervals.

x = -M/2*T0:0.01:M/2*T0; 

syms t n y a
y=sym('Heaviside(t+a)')*A*2-sym('Heaviside(t-a)')*A*2;
y=subs(y,a,tau/2);
y=simple(y);

%plot the input square wave
figure1 = figure(1);
axes1 = axes('FontSize',14,'Parent',figure1);
box(axes1,'on');
hold(axes1,'all');
ylim(axes1,[dc-2*A dc+2*A]);
grid;
title(['Square wave']);
xlabel('t (seconds)'); ylabel('Amplitude');

xx = (x>=0).*(x-ts-fix((x-ts+T0/2)/T0).*T0) + (x<0).*(x-ts-fix((x-ts-T0/2)/T0).*T0);
yy = double(subs(y,t,xx))-A+dc;
yy(isnan(yy)) = 0;
plot(x,yy,'LineWidth',2,'color','b'); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% closed-form of Cn (KEY EQUATION) %%%%%%%%%%
%Compute the closed-form of Cn; i.e.,
%C0 = 1/T0 * integral(f(t)) over [-T0/2, T0/2]
%Cn = 1/T0 * integral(f(t)*exp(-j*n*w0*t)) over [-T0/2, T0/2]
C0=int(y,t,-T0/2,T0/2)/T0;
Cs=int(y*exp(-j*w0*n*t)/T0,t,-T0/2,T0/2);
C0=simple(C0);
Cs=simple(Cs);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% closed-form of Cn (KEY EQUATION) %%%%%%%%%%

%Figure 2, plot amplitude spectrum
figure2 = figure(2);
axes2 = axes('FontSize',14,'Parent',figure2);
box(axes2,'on');
hold(axes2,'all');
title(['Amplitude Spectrum']);
xlabel('n'); ylabel('Amplitude');

%Compute Cn and plot its amplitude
for k=-Nf:1:Nf
    if k==0
        c = abs(double(C0));             %C0
    else
        c = double(subs(Cs,n,k));   %Cn, substitue n=k into the above closed-form of Cn
    end

    stem(k, abs(c),'color','b','MarkerSize',5);                  %plot
    
end

%Figure 3, plot phase
figure3 = figure(3);
axes3 = axes('FontSize',14,'Parent',figure3);
box(axes3,'on');
hold(axes3,'all');
title(['Phase Spectrum']);
xlabel('n'); ylabel('Phase (degrees)');

%Compute Cn and plot its angle
for k=-Nf:1:Nf
    if k==0
        c = abs(double(C0));             %C0
    else
        c = double(subs(Cs,n,k));   %Cn, substitue n=k into the above closed-form of Cn
    end

    stem(k, sign(k)*angle(c)*180/pi,'color','b','MarkerSize',5);    %plot
    
end