% Recontruct step-by-step, by user prompt, the waveform of some Periodic 
% functions from their Fourier Series coefficients
%
%Enter the value of the string variable "waveform" below as follows:
%sq--for a square wave
%6p--for a 6-pulse square wave
%rs--for a rectified sinewave
%sw--for a sawtooth
%tr--for a triangular wave
%
%Hit left mouse button to see the stages of reconstruction; Rght mouse to
%stop

close all; clc; clear all

waveform='sq';

%Define time
t=-1.5:0.005:1.5;%Plot time points
tm=mod(t+1/2,1)-1/2;%normalized time points to the intervale [-1/2 1/2]

switch lower(waveform)
    case 'sq'
        f=-1*(((-0.5<=tm)&(tm<=-0.25))|((0.25<=tm)&(tm<=0.5)))...
            +1*((-0.25<=tm)&(tm<=0.25));
        r=[0:1:7];
        n=2*r+1;%Harmonic number
        A=(4/pi)*(1./(((-1).^(3*r+2)).*n));%Fourier cosine term amplitudes
        B(r+1)=0;%Fourier sine term amplitudes
    case '6p'
        f=-1*( ((-0.5<=tm)&(tm<=-1/3))|((1/3<=tm)&(tm<=0.5)) )+1*((-1/6<=tm)&(tm<=1/6));
        r=0:1:5;
        n=[6*r+1 6*r(2:end)-1];n=sort(n);%Harmonic order
        A=((n-6*round(n/6))).*(2*sqrt(3)/pi).*(1./n);%Fourier cosine term amplitude
        B(2*r+1)=0;%Fourier sine term amplitude
    case 'rs'
        f=abs(sin(2*pi*t));
        r=[0:1:7];
        n=2*r;%Harmonic number
        A=-(4/pi)*1./((n-1).*(n+1));%%Fourier cosine term amplitudes
        A(1)=A(1)/2;%Average value
        B(r+1)=0;%Fourier sine term amplitudes
    case 'sw'
        f=2*tm;
        r=0:1:7;
        n=r+1;%Harmonic number
        A(r+1)=0;%Fourier term cosine amplitude
        B=((-1).^r).*(2/pi)./n;%Fourier sine term amplitude
    case 'tr'
        f=((-0.5<=tm)&(tm<=0)).*(tm+0.25)-1*((0<tm)&(tm<=0.5)).*(tm-0.25)
        r=0:1:7;
        n=2*r+1;%Harmonic Number
        A=(2/pi/pi)./(n.^2);%Fourier cosine term amplutide
        B(r+1)=0;%Fourier sine term amplitude
end

figure %Plot time waveform
b=0;k=0;F=0;
tmn=min(t);tmx=max(t);fmn=1.3*min(f);fmx=1.3*max(f);
Abmx=max(abs(fmn),abs(fmx));

subplot(2,1,1)
plot(t,f,'k','LineWidth',1.5)
line([tmn tmx],[0 0], 'Color','k','LineStyle','--','LineWidth', 1.2)
line([0 0],[fmn fmx],'Color','k','LineStyle','--','LineWidth', 1.2)
axis([tmn tmx fmn fmx])
title('Waveform v. Time')

subplot(2,1,2)
plot(t,f-F,'b','LineWidth',1.5)
line([tmn tmx],[0 0], 'Color','k','LineStyle','--','LineWidth', 1.2)
line([0 0],[-Abmx Abmx],'Color','k','LineStyle','--','LineWidth', 1.2)
axis([tmn tmx -Abmx Abmx])
title('Error v. Time')
xlabel('Time, normalized to 1 period')    

[x,y,b]=ginput(1);
while b==1
    k=k+1;
    F=F+A(k)*cos(2*pi*n(k)*t)+B(k)*sin(2*pi*n(k)*t);%Harmonic term accumulation 
    subplot(2,1,1)
    plot(t,f,'k',t,F,'r','LineWidth',1.5)
    line([tmn tmx],[0 0], 'Color','k','LineStyle','--','LineWidth', 1.2)
    line([0 0],[fmn fmx],'Color','k','LineStyle','--','LineWidth', 1.2)
    axis([tmn tmx fmn fmx])
    title(['Waveform v. Time: adding n=' num2str(n(k)) ' harmonic order'])
    
    subplot(2,1,2)
    plot(t,f-F,'b','LineWidth',1.5)
    line([tmn tmx],[0 0], 'Color','k','LineStyle','--','LineWidth', 1.2)
    line([0 0],[-Abmx Abmx],'Color','k','LineStyle','--','LineWidth', 1.2)
    axis([tmn tmx -Abmx Abmx])
    title('Error v. Time')
    xlabel('Time, normalized to 1 period')
    [x,y,b]=ginput(1);
    if k==max(r)+1; b=0;end
end

figure %Harmonic Amplitude plot
subplot(2,1,1)
bar(n,sqrt(A.^2+B.^2),0.25,'r')
title(['Harmonic Amplitude v. Harmonic Number for the ', waveform, ' waveform'])
ylabel('Harmonic Amplitude, Normalized')

subplot(2,1,2)
bar(n,atan(-B./A),0.25,'r')
title(['Harmonic Phase v. Harmonic Number for the ', waveform, ' waveform'])
xlabel('Harmonic Number')
ylabel('Harmonic Phase, rad')