% ### Swave3.m ###   03.02.10

% plots a user-specified # of terms for the Fourier series expansion of a
% square wave; also animates the plot

clear
clf

% --------------

N= 10;   % number of non-zero terms (including the DC term; must have N >= 0)

% --------------


% define time
t= linspace(-3.5,3.5,1000);

% define square wave (brute force!)
for nn=1:size(t,2)
    if t(nn)>=-3&&t(nn)<-2 || t(nn)>=-1&&t(nn)<-0 || t(nn)>=1&&t(nn)<2 || t(nn)>=3&&t(nn)<4
        f(nn)= 0;

    else
        f(nn)= 1;
    end

end

plot(t,f,'o')
axis([-3.5 3.5 -0.25 1.25])
grid on; hold on;
xlabel('t')
ylabel('f(t)')



% now, compute the Fourier series with the specified # of terms

FS= 0;  % initially set array containing Fourier series values

for nn=1:N+1

    % solution below stems frm ch.10.5 in Hughes-Hallet
    if nn==1
        FS= 1/2*ones(1,size(t,2));    % DC term
        
        % 'animate' things getting plotted
        clf
        plot(t,f,'o')
        axis([-3.5 3.5 -0.25 1.35])
        grid on; hold on;
        xlabel('t')
        ylabel('f(t)')
        plot(t,FS,'r-','LineWidth',2)
        legend('Square wave',['Fourier Approximation',' (',num2str(nn-1),'/',num2str(N),')']);
        pause(0.5) 
        
    else
        cnt= nn-2;  % dummy index
        FS = FS + (2/((2*cnt+1)*pi))*sin(pi*(2*cnt+1)*t); % n'th Fourier term

        % 'animate' things getting plotted
        clf
        plot(t,f,'o')
        axis([-3.5 3.5 -0.25 1.35])
        grid on; hold on;
        xlabel('t')
        ylabel('f(t)')
        plot(t,FS,'r-','LineWidth',2)
        legend('Square wave',['Fourier Approximation',' (',num2str(nn-1),'/',num2str(N),')']);
        pause(0.5)
        
    end

end

legend('Square wave',['Fourier Approximation',' (',num2str(N),' non-zero terms included)'])