% FUNCTION: Harmon(y,t,w)
%
% Usage:
%     [Mag, Phase, Mean ] = harmon(Y,T,W)
%
% Input:
%
%	Y	Dependent variable containing N points.
%	T	Independent variable containing N points.
%	W	Array containing M distinct frequencies.
%
% Output:
%
%	Mag	Magnitude of each harmonic constituent.
%	Phase	The phase (in radians) of each constituent.
%	Mean	An estimate of the mean value.
%
% Theory:
%
% This function performs a least-squares harmonic
% fit to irregularly-spaced data.  We solve the
% least-squares problem assuming that the values
% Y(i) can be modeled as a sum of sines and cosines
% along with a single constant term:
%
%     f(t) = A_0 + \sum_{n=1}^M {
%                A_n \cos(\omega t) + B_n \sin(\omega t)
%                }
%
% Once f(t) has been defined this way, it is fairly
% straightforward to derive the 2M+1 equations
% necessary to solve for all the coefficients.
% This exercise is left to the reader...

function [Mag, Phase, Mean ] = harmon(y,t,w)

%%%%%%
% Making sure that I'm getting column vectors
% and not row-vectors...

Y = y ;
if  size(Y,1) == 1 
    Y = Y' ;
end

T = t ;
if  size(T,1) == 1 
    T = T' ;
end

%%%%%%
% 'M' is the total number of frequencies...
% 'N' is the total number of data points...

M = length(w) ;
N = length(Y) ;

%%%%%%
% I need to solve a matrix system of the form
%
%    A * X = B
%
% This is the most compact way I found to build
% the 'A' and 'B' matrices...

A  = ones(N,1) ;

for i = 1:M
    A = [ A cos( w(i)*T) sin( w(i)*T) ] ;
end

B = A' * Y ;

A = A' * A ;

%%%%%%
% This solves the system of equations.
% It's rather anti-climatic if you ask me...

X = A \ B ;

Mean = X(1) ;

%%%%%%
% It makes more sense to express things as
%     A_n \cos( \omega t + \phi_n)
% compared to
%     A_n \cos( \omega t )  + B_n \sin( \omega t )

k = 2:2:2*M ;

Mag   = sqrt( X(k).^2 + X(k+1).^2 ) ;
Phase = atan2( X(k+1), X(k) ) ;

clear A B k M N T X Y

