clear all
close all
% Fourier analysis for a square wave signal

l_w=2;
figure
hold on
t=linspace(0,2,500);
y2=[ones(125,1);zeros(125,1);ones(125,1);zeros(125,1)];


plot(t,y2,'--g','linewidth',2);
ylim([0 1.1])
grid on 

figure
hold on
t=linspace(0,2,500);
y2=[ones(125,1);zeros(125,1);ones(125,1);zeros(125,1)];
% Single harmonic approximation
y=1/2+2/pi*sin(2*pi*t);
subplot(3,1,1)
plot(t,y,t,y2,'--r','linewidth',l_w);
title('Single harmonic approximation','FontSize',12)

% Two harmonics approximation
y=1/2+2/pi*sin(2*pi*t)+2/pi/3*sin(3*2*pi*t);
subplot(3,1,2)
plot(t,y,t,y2,'--r','linewidth',l_w);
title('Two harmonics approximation','FontSize',12)
ylabel('Amplitude','FontSize',16)

% Five harmonics approximation
y=1/2+2/pi*sin(2*pi*t)+2/pi/3*sin(3*2*pi*t)+2/pi/5*sin(5*2*pi*t)+2/pi/7*sin(7*2*pi*t)+2/pi/9*sin(9*2*pi*t);
subplot(3,1,3);
plot(t,y,t,y2,'--r','linewidth',l_w)
title('Five harmonics approximation','FontSize',12)
xlabel('time [s]','FontSize',16)


% Frequency domain analysis of the effect of a square wave input on a low
% pass filter
s=tf('s');
Gs=1/(s^2+s+1);
w=logspace(-1,2,500);
[M,P]=bode(Gs,w);
figure
subplot(2,1,1);
semilogx(w,20*log10(squeeze(M)),'LineWidth',2);
axis([0.1,10,-40,10])
grid on
title('Magnitude frequency response of the system','FontSize',12)
ylabel('Amplitude [dB]','FontSize',14)

% Plot the frequency spectrum of the signal
subplot(2,1,2)
wv=[1,3,5,7,9,11];
Mv=[2/pi,2/(3*pi),2/(5*pi),2/(7*pi),2/(9*pi),2/(11*pi)]; % 2/(n*pi)
Mv_dB=20*log10(Mv);
semilogx(wv,Mv_dB,'o','LineWidth',5,'MarkerSize',5);
axis([0.1,10,-40,10])
hold on
grid on
for k=1:length(wv)
    plot([wv(k),wv(k)],[Mv_dB(k)-40,Mv_dB(k)],'LineWidth',2);
end
ylabel('Amplitude [dB]','FontSize',14)
title('Magnitude of the first five Fourier coefficients of a square wave signal with T=2\pi','FontSize',12)
xlabel('\omega [rad/s]','FontSize',14)


% Output of the system against pulse wave and two-terms approximation from
% the Fourier expansion
% Pulse wave
[uv1,tv1]=gensig('square',2*pi,50,5/1000);
uv1=1-uv1;
% Average value plus first harmonic
[uv2,tv2]=gensig('sin',2*pi,50,5/1000);
uv2=2/pi*uv2+1/2;


yv1=lsim(Gs,uv1,tv1);
yv2=lsim(Gs,uv2,tv2);

figure
hold on
grid on
plot(tv1,yv1,'--r',tv2,yv2,'-','LineWidth',2)
legend('square wave output','first harmonic approx.','Orientation','horizontal');
ylabel('Output','FontSize',14)
xlabel('time [s]','FontSize',14)



tu=linspace(0,4*2*pi,4*1000);
u2a=[ones(500,1);zeros(500,1)];
u_2=[u2a; u2a; u2a; u2a;];
% Single harmonic approximation
u_1=1/2+2/pi*sin(tu);
yu1=lsim(Gs,u_1,tu);
yu2=lsim(Gs,u_2,tu);

figure
hold on
grid on
plot(tu,yu1,'-',tu,yu2,'--r','LineWidth',2)
legend('first harmonic approx.','square wave output','Orientation','horizontal');
ylabel('Output','FontSize',14)
xlabel('time [s]','FontSize',14)


figure % plot the two ways to define a square wave and the correspondig approx.
subplot(2,1,1)
plot(tv1,uv1,'--r',tv1,uv2,'-');
legend('square wave','first harmonic approx.','Orientation','horizontal');
subplot(2,1,2)
plot(tu,u_2,'--r',tu,u_1,'-');
legend('square wave','first harmonic approx.','Orientation','horizontal');