%FourierApprox inputs a function handle f on the interval [-pi,pi] which
%is extended periodically to [-3pi,3pi]. It then approximates this periodic
%function by Fourier polynomials of order M,2M,...floor(N/M)M. The results
%are plotted and each sucessive approximation is diplayed by "Press any
%key".


%By Evan Randles for the I.H.S. math seminar

function[]=FourierApprox(f,N,M)
x=[-pi:.01:pi];                             %Interval on which f is defined
fvec=f(x);                                  %Discretization of f
z=[-3*pi:.01:3*pi];                         %f is extended periodically
g=[fvec(1:1:length(fvec)-1),fvec(1:1:length(fvec)-1),fvec]; 
plot(z,g)                                   %plot of periodic extension
legend('The original function')
pause
hold on
for i=[1:N]                                 %Computation of Fourier Coefs.
    b(i)=(pi)^(-1)*quad(@(x)SinInt(x,i),-pi,pi);   %Uses quad function
    a(i)=(pi)^(-1)*quad(@(x)CosInt(x,i),-pi,pi);   
end
legendinfo=cell(1,floor(N/M)+1);            %Initialization of string for  
legendinfo{1}='The original function';      %use in legend
for i=2:floor(N/M)+1                        
    legendinfo{i}='';
end
legendindex=2;
yDummy=(2*pi)^(-1)*quad(f,-pi,pi)*ones(1,length(z));
for k=[1:N]                                %Construction of Fourier Poly
        y(k,:)=yDummy(:)+a(k)*cos((k)*z(:))+b(k)*sin(k*z(:));
        yDummy=y(k,:);
        if mod(k,M)==0;
            plot(z,y(k,:),'Color',rand(1,3)) %Plot of Fourier Poly w/legend
            legendinfo{legendindex}=['Fourier Polynomial of degree ' num2str(k)];
            legend(legendinfo(1:legendindex)) %Updating Legend
            legendindex=legendindex+1;
            pause
        end 
         
end

hold off
    function y=CosInt(x,n)   %Cosine integrand for computation of 
        y=f(x).*cos(n*x);    %Fourier coefficients for quad function
    end
        function y=SinInt(x,n)  %Sine integrand for computation of 
        y=f(x).*sin(n*x);       %Fourier coefficients for quad function
    end
end