function [] = fseriesdemo(f,b,NN)
% FSERIESDEMO Plots a function f and it's Fourier series approximation.
% FSERIESDEMO(f,b,N) plots the function f and it's N term Fourier series
% approximation from zero to b.  FSERIESDEMO uses the FFT to approximate
% the Trigonometric Fourier Series coefficients.  FSERIESDEMO only works
% for the function from 0 to b.  Input N must be an integer >=2.
% Example:
%
%           f = @(x) sin(x)./x+(cos(x)).^2-exp(x/8)+x.^(x/20);
%           fseriesdemo(f,20,100)  % Approx. f on (0,20) with 100 terms.
%
%
% Author:  Matt Fig
% Contact:  popkenai@yahoo.com

% First some argument checking.
if nargin<3
   error ('Not enough arguments. See help.')
end

if ~isreal(NN) || ~isreal(b) || b<=0 || NN < 2 || floor(NN)~=NN
    error ('2nd and 3rd args must be positive real scalars. See help.')
end

n = 2*NN;
T = b-eps+(b-eps)/(n-2); % Get the appropriate endpoint for fft.
x = linspace(eps,T,n+1); % Use this x in fft.
x(end) = [];   % But not the endpoint.
fun = f(x);  % Get discrete points.
FUNC = fft(fun);  % Call fft
FUNC = [conj(FUNC(NN+1)) FUNC(NN+2:end) FUNC(1:NN+1)]/n; % Shift/scale.
A0 = FUNC(NN+1);   % First cosine term.
AN = 2*real(FUNC(NN+2:end));  % General cosine term.
BN = -2*imag(FUNC(NN+2:end));  % General sine term.


%My Version
function [a0 ak bk]=fcoeffs_fft(profile,angles,Ncoeffs)
%FCOEFFS_FFT Fourier coefficients (trigonometric) via Fast Fourier Transform (FFT)
%  [a0,ak,bk]=fcoeffs_fft(profile,angles,Ncoeffs) determines the Fourier coefficients of
%  irregulary spaced angles and profile.
%  Ncoeffs is the number of Fourier coefficients to return
%  Oversamples input via a (periodic) cubic spline interpolation
%  Returns DC component a0, cosine terms ak and sine terms bk


%Oversampling the data points so the total number of points is Ntheta
Ntheta = 1024;
theta = linspace(-pi,pi,Ntheta);
delta_theta = 2*pi/Ntheta;

%This defines a temporary vector that "wraps around" pi to -pi
theta_over = -pi:delta_theta:(pi + 5*delta_theta);

%Creating the function to which the fft can be applied
pp = spline(angles,profile);
zeta_temp= ppval(pp,theta_over);
%Overwrite here
zeta = zeros(size(theta));
zeta(1:Ntheta) = zeta_temp(1:Ntheta);
zeta(1:6)=zeta_temp(Ntheta:Ntheta +5);

zetahat=fft(zeta);

a0 = zetahat(1)/Ntheta;                                      %DC component

for ik = 1:Ncoeffs
    ij = ik +1;       %index for fft

ak(ik)=2*real(zetahat(ij))/(Ntheta);                           %cosine terms
bk(ik)=-2*imag(zetahat(ij))/(-1.*Ntheta);                       %sine terms
end