%**************************************************************************
%Initialization and Setup
clear freq x w t n ...  %reset all variables to zero
        i f numharms timepts harmonics Vdc HAtotal HA A        

%Harmonic Calculation Variables
t=0:5e-6:1e-3;              % Time; this is the "x" value in the plot3 command.
A=2;                        % Amplitude
f=1000;                     % Frequency
w=(2*pi)*f;                 % Omega
Vdc=0;                      % DC component
n=11;                       % Last harmonic number

%x, y, and z-plane Matrix setup and plot variables
timepoints=length(t);           % Number of data points in each harmonic plot
harmonics=(1:2:n);              % Harmonic numbers
numharms=length(harmonics);     % Number of harmonics in frequency domain
                                %   numharms != harmonics because a square
                                %   wave only uses odd harmonics in the
                                %   Fourier series
il=1;                           % Index location; used for harmonic array HA
freq=ones(numharms,timepoints); % This is an array of row vectors.  All 
                                %   elements in each row contain the same
                                %   value, which corresponds to the
                                %   harmonic number for that row.  These
                                %   are the "y" values in the plot3
                                %   command.
HA=zeros(numharms,timepoints);  % Harmonic Amplitude; this is a 2-D array
                                %   containing the voltage magnitude for
                                %   each point in time (one dimension) for
                                %   each harmonic (2nd dimension).  These
                                %   are the "z" values in the plot3
                                %   command.

%**************************************************************************
% Initialize each row vector with its corresponding harmonic number.
for i=(1:numharms)
    freq(i,:)=harmonics(i);
end
%**************************************************************************
% Initialize getframe and movie variables; create avi output file.
num=1;
nframes = 1;
Frames = moviein(nframes);      % Movie frames
mov = avifile('D:\Data\Research\StudentProjects\FourierConceptIllustration\TimeFreqRotation.avi','compression','Cinepak');
%**************************************************************************
% Calculate amplitude value for each time point in each harmonic
for x=1:2:n                             %will only calculate odd harmonics
    HA(il,:)=((2*A)/(pi*x))*sin(x*w*t); %Harmonic Amplitude calculation
    il=il+1;
end
%**************************************************************************
% Create 3-D plot of Fourier series
cla;
il=1;
for x=1:2:n
    hold on
    %plot of harmonic series
    plot3(t,freq(il,:),HA(il,:),'LineWidth',2)  
    il=il+1;              
end
    xlabel('Time (s)'), ylabel('Frequency (kHz)'), zlabel('Voltage (V)')
    xlim([-0.05e-3 1.05e-3]);
    ylim([-1 n+1]);
    zlim([-1.5 1.5]);
%**************************************************************************
%View rotation commands
slowdown=3; %Take this many extra frames per view instead of just one early 
            %   in rotation to slow down video.
for c=90:-2:0
    view([c 0]);
    if c>74     % Take extra frames until this rotation to slow down the video.
        for localcnt=1:1:slowdown
            Frames(num)=getframe(gcf);
            num=num+1;
            pause(0.2); % The pause statements are for effect when plotting 
                        %   from within MATLAB.
        end
    end
    pause(0.2);         % For effect when plotting from within MATLAB.
    Frames(num)=getframe(gcf);              
    num=num+1;
end
for localcnt=1:1:20     % Take extra frames to "pause" video looking at 
                        %   overlaid harmonics.
    Frames(num)=getframe(gcf);
    num=num+1;
end
pause(2);               % For effect during MATLAB plot.
%**************************************************************************
% Change to square wave
xsquare =[0; 0; 0.5e-3; 0.5e-3; 1e-3; 1e-3];
ysquare =[0; 1;   1;      -1;    -1;    0];
cla reset
plot (xsquare,ysquare,'LineWidth',2)
xlabel('Time (s)'), ylabel('Voltage (V)')
xlim([-0.05e-3 1.05e-3]);
ylim([-1.5 1.5]);
for localcnt=1:1:20  % Take extra frames to "pause" video looking at square wave result.
    Frames(num)=getframe(gcf);
    num=num+1;
end
%**************************************************************************
% Finish movie and close avi file.
mov = addframe(mov,Frames);
mov = close(mov);
%**************************************************************************