function [out, dc_term] = fourier_analysis(signal, num_harmonics)
%function [out, dc_term] = fourier_analysis(signal, num_harmonics)
%
% FOURIER_ANALYSIS: Performs Fourier series analysis on one period
%                   of a real (non-complex) signal.
%
% Input Parameters:
%    signal        -- One period of a periodic signal to be analyzed
%    num_harmonics -- The number of positive-frequency coefficients to
%                     calcuate.  This parameter must be less than or 
%                     equal to one-half the length of signal. (Default: 10)
%
% Output Parameters:
%    out     -- the coefficients resulting from the Fourier series analysis.
%               Out contains the complex amplitudes for a cosine reconstruction
%               of the signal from 1 to num_harmonics.  
%    dc_term -- the DC offset term (or "zeroth Fourier coefficient")
%

roundoff_thresh = 1e-14;

% Check input parameters
if ~exist('num_harmonics','var') | isempty(num_harmonics)
    num_harmonics = 10;
end

if size(signal,2) > size(signal,1) 
    signal = signal';
end

if num_harmonics > length(signal)/2
    error(['Cannot extract more than ' num2str(floor(length(signal)/2)) ' harmonics.']);
end

% Set up some useful arrays
N = length(signal);
n = (0:(N-1));
out = zeros(1,num_harmonics);

dc_term = sum(signal)/N;
for k = 1:num_harmonics;
    ce = exp(j*k*n/N*2*pi)';            % Build the CE signal
    out(k) = sum(ce.*signal)/N;         % Do the correlation
end

% Because we correlated with complex exponential signals, we're
% really off by a factor of 2 from the required cosine amplitude
% (since a cosine is the sum of two CE signals). 
%  -- Note that this follows the conventions of "DSP First!", ch3
out = out*2;

% The component at the Nyquist frequency (fs/2) has no
%  equivalent negative frequeny.  Since we're representing all
%  other signals as pairs, the Nyquist frequency needs to have
%  its amplitude cut in half.
if num_harmonics == length(signal)/2
    out(:,end) = out(:,end)/2;
end

% Let's get rid of coefficients that are less than roundoff_thresh
r = real(out);
r(find(abs(r) < roundoff_thresh)) = 0;
im = imag(out); 
im(find(abs(im) < roundoff_thresh)) = 0;

out = r + i*im;
