%% Fourier Series
%
% In this lab, we will learn how to create, manipulate, and plot Fourier
% series. 
%
%% Reference
%
% This is Computer Lab 4.2 *fourierseries* from the text
% <http://www.stthomas.edu/wavelets Discrete Wavelet Transformations: An
% Elementary Approach with Applications> by Patrick J. Van Fleet, Wiley
% 2008.

%% Matlab Commands
%
% The matlab commands used in this lab are: |mkpp|, |ppval|, |trapz|, and 
% |conv|.
%
%% DiscreteWavelets Commands
%
% This lab makes use of the 
% <http://www.stthomas.edu/wavelets DiscreteWavelets Toolbox>.  
% Commands used in this lab are: |FiniteFourier|, |Cs|, |Sn|.
%
%% Due Date
% 
% This lab is due on *DUE DATE*.  You may work with a partner if you wish
% to do so.
%
%% Instructions
%
% Place this m-file in a folder with all other lab m-files and save it.
% Access the m-file and then click on File -> Publish to HTML from the menu
% above.  This will produce a nice html file that can be easily viewed with
% the Matlab web browser.  
%
% Edit this m-file using the Matlab editor and insert your answers to the
% problems below.
%
%% Computing Fourier Series in Matlab
%
% Let's look at an example to understand how to work with Fourier series in
% Matlab
%
%% 
% We consider the 2pi-periodic function f(w) where f(w) = w on the interval
% (-pi, pi).
%
% Here are three periods of f(w):
%
polycoeffs=[1 -pi];
plotrange=-3*pi:0.1:3*pi;
piecewise=mkpp([-3*pi -pi pi 3*pi],[polycoeffs;polycoeffs;polycoeffs]);
plot(plotrange,ppval(piecewise,plotrange),'k-')
%% Computing the Fourier Coefficients
%
% Since f(w) is an odd function, we know the Fourier coefficients must be
% imaginary and from Example 4.4 on Page 113, we have a formula for 
% computing the Fourier coefficients c_k.  In Matlab we can build the
% following coefficient vector numerically using the trapezoidal rule.
%
% First set the number of points
%
np=5;
%
% Now create the sample points and a vector to hold the coefficients
%
x = linspace(0,pi,np);
c = zeros(1,np);
%
% Now we can create the Fourier coefficients
%
for k=1:np
    y=x.*sin(k*x);
    c(k)=-i*trapz(x,y)/pi;
end
%% Create and Plot the Series
%
% Now we can build the sine series and plot it.  First we must modify the
% Fourier coefficients as done in Example 4.4.
%
d=2*i*c;
%
% Now we create the partial series - we make it for three periods
%
fseries=0;
for k=1:np
    fseries=fseries+d(k)*sin(k*plotrange);
end

plot(plotrange,[ppval(piecewise,plotrange); fseries]);
set(gca,'XTick',-3*pi:pi:3*pi)
set(gca,'XTickLabel',{'-3*pi','-2*pi','-pi','0','pi','2*pi','3*p'})
title(['Number of coefficients: ' int2str(np)])
%% Finite Length Fourier Series
%
% We will often have occasion to start with some finite list of terms
% a_0,...,a_n say, and construct a Fourier series from then.  
% DiscreteWavelets has a useful command called |FiniteFourier| for 
% constructing such a series.  The module takes a vector of sample points,
% the list of coefficients, and a starting index.
%
h=[1/4 0 1/2 0 1/4]; % Fourier coefficients for the function cos^2(w/2)
w=-pi:.001:pi; %Create sample points
plot(w,FiniteFourier(w,h,-2)); %We need to start the series at k=-2.
set(gca,'XTick',-pi:pi/2:pi)
set(gca,'XTickLabel',{'-pi','-pi/2','0','pi/2','pi'})
title(['Fourier series for cos^2(w/2)'])
%%
% Some Fourier series make sense to plot in the complex plane, but by and
% large, we will plot only the moduli of Fourier series.
%
v=[2 1];
w=-pi:.001:pi;
%Giving plot one complex-valued argument results in a parametric plot.
plot(FiniteFourier(w,v,0)); 
title('The Fourier series H(w) = 2+e^{iw}')
set(gca,'DataAspectRatio',[1 1 1])
%% Fourier Series from Sin and Cos
%
% Since cos(w) and sin(w) can be written in terms of complex exponentials,
% we can often expand them as such and form Fourier series.  
% DiscreteWavelets contains two functions, |Cs| and |Sn|, that give the
% Fourier coefficients for sin(w) and cos(w) as vectors.
%
Cs()
%
% The output here is a three vector - you need to remember that the
% elements correspond to e^{-iw}, 1, and e^{iw}, respectively.
%
% We can make a Fourier series for cos(w) and plot it.  Note that the
% starting index for |FiniteFourier| is -1.
%
h=Cs();
w=-pi:.001:pi;
plot(w,FiniteFourier(w,h,-1));
set(gca,'XTick',-pi:pi/2:pi)
set(gca,'XTickLabel',{'-pi','-pi/2','0','pi/2','pi'})
title('The Fourier series cos(w)')
%
%%
% We can create and plot sin(w) in the same manner:
%
h=Sn();
w=-pi:.001:pi;
plot(w,FiniteFourier(w,h,-1));
set(gca,'XTick',-pi:pi/2:pi)
set(gca,'XTickLabel',{'-pi','-pi/2','0','pi/2','pi'})
title('The Fourier series sin(w)')
%
%%
% What is the Fourier series for  sin^3(w)?  Plot the modulus of this
% series.
%
% We use the Matlab command |conv| to solve this problem:
h=conv(Sn(),Sn());
h=conv(Sn(),h)
%
% We know the indices of this series run from -3 to 3.  So we create and
% plot the modulus of the Fourier series accordingly.
w=-pi:.001:pi;
plot(w,abs(FiniteFourier(w,h,-3)));
set(gca,'XTick',-pi:pi/2:pi)
set(gca,'XTickLabel',{'-pi','-pi/2','0','pi/2','pi'})
title('The modulus of the Fourier series sin^3(w)')

%% Exercises
%
% 
% *Problem 1* 
%
% Suppose g(w) = pi^2 - w^2 for -pi <= w <= pi$$ and 0 otherwise.  Now form 
% the 2 pi - periodic function f(w) = \sum_k g(w-2*k*pi).  A plot of three
% periods of f(w) appears below:
%
polycoeffs=[-1 2*pi 0];
plotrange=-3*pi:0.1:3*pi;
piecewise=mkpp([-3*pi -pi pi 3*pi],[polycoeffs;polycoeffs;polycoeffs]);
plot(plotrange,ppval(piecewise,plotrange),'k-')
set(gca,'XTick',-3*pi:pi:3*pi)
set(gca,'XTickLabel',{'-3*pi','-2*pi','-pi','pi','2*pi','3*pi'})
%
% Use Matlab to find the Fourier coefficients for f(w).
%
%%
% Place your answer here.  To type an answer, start a line with %.  Matlab
% commands are entered as usual.

%%
% *Problem 2* 
%
% Find the Fourier coefficients of the function m(w) = f ' (w) where f(w)
% is given in Problem 1.  Hint: There is a much faster way for finding 
% these coefficients than direct integration!! 
%
%%
% Place your answer here.  To type an answer, start a line with %.  Matlab
% commands are entered as usual.
%
%%  
% Once you have the Fourier coefficients, use them to construct a  series
% in terms of sin(kw) that we can use to plot partial sums that approximate 
% m(w).
%
%%
% Place your answer here.  To type an answer, start a line with %.  Matlab
% commands are entered as usual.
%

%%
% *Problem 3* 
%
% Find the coefficients (real-valued, non-zero) for a length 4 Fourier
% series H(w) so that H(0)=1 and H(pi)=0.  There is no restriction on the 
% indices of the four coefficients.  They could be 0, 1, 2, 3,  or 
% -2, 5, 9, 15, etc.  - your choice.
%
%%
% Place your answer here.  To type an answer, start a line with %.  Matlab
% commands are entered as usual.

%% 
% *Problem 4*
%
% Find the Fourier series coefficients for the 2pi - periodic function
% f(w) = cos^3(w)*sin^2(2w).
%%
% Place your answer here.  To type an answer, start a line with %.  Matlab
% commands are entered as usual.

%%
% *Problem 5 (Challenge)*
%
% Find the coefficients (real-valued) for a finite length Fourier series 
% H(w) such that H(0)=1, and H'(0) = H(pi) = H'(pi) = 0.  What is the 
% minimum number of nonzero terms that allows a unique solution?  
% What is this solution?
%
% Hint:  The conditions above lead to a linear system.   Find that system
% by hand and then use Matlab to solve it.  You will have to experiment
% with the number of terms in H(w) to find the right length that yields a
% unique solution.  
%%
% Place your answer here.  To type an answer, start a line with %.  Matlab
% commands are entered as usual.


%%
close all;
displayEndOfDemoMessage(mfilename)