% Square Wave from Sine Waves
% The Fourier series expansion for a square-wave is made up of a sum of odd
% harmonics.  
%
%   Modified by Kevin D. Donohue (donohue@engr.uky.edu) August 29, 2011
%    adapted from Matlab's xfourier.m script.

nhar = 15; %  Number of harmonics in expansion (only odd harmonics are used)
fs = 1024;  % sampling freqency in Hz
tinc = 1/fs;  % sampling increment
t = 0:tinc:1;   %  Create time sample points over a 1 second interval
y = zeros(nhar,length(t)); % initalize matrix to save cumulative reconstructions 
x = zeros(size(t));  %  Initalize accumulator for waveform
%  Loop to add up odd harmonics
colr = ['b', 'r', 'g', 'k', 'c'];  % Colors for multiple waveforms on same plot 
for kl=1:nhar
    k = 2*(kl-1)+1;  % convert to odd number
    figure(2)  % figure to plot individual sine waves
    set(gcf,'Position', [567   429   437   247])
    sad = (4/(k*pi))*sin(2*pi*k*t);  %  Individual harmonic
    plot(t,sad,colr(mod(kl-1,length(colr))+1)); hold on
    xlabel('Seconds')
    ylabel('Amplitude')
    title('Individual Components')
    x = x + sad;  %  Add up harmonics
    y(kl,:) = x;  % Save culumative waveform at this stage
    figure(1)  %  Figure to plot cumulative waveform
    set(gcf,'Position', [62   429   437   247])
    plot(t,x)
    xlabel('Seconds')
    ylabel('Amplitude')
    title('Composite Waveforms')
  %  pause   %  Pause to update graphics and look at plots
end
figure(2)
hold off

%pause

% Here is a 3-D surface representing the gradual transformation of a sine wave
% into a square wave.

surf(y);
shading interp
axis off ij
title('The building of a square wave: Gibbs'' effect')