% This function is part of the NMSM Pipeline, see file for full license.
%
% Computes up to 7th degree polynomial plus Fourier fitted values given
% polynomial and Fourier coefficients saved in a column vector, as found by
% polyFourierCoefs.m. This function is intended to by vectorized.
%
% (Array of double, double, Array of double, double, logical)
% -> (Array of double)
% Computes polynomial and Fourier fitted values from coefficients.
function fit = polyFourierCurves(coefs,f,t,degree,derivs)
% Calculate frequency in radians/second
w = 2*pi*f;
% Calculate cosine and sine harmonic arrays
ncoefs = size(coefs,1);
nharmonics = (ncoefs-(degree+1))/2;
npts = size(t,1);
% Calculate y
a = zeros(npts,ncoefs);
switch degree
case 0
a(:,1) = 1.0;
case 1
a(:,1) = 1.0;
a(:,2) = t;
case 2
t2 = t.*t;
a(:,1) = 1.0;
a(:,2) = t;
a(:,3) = t2;
case 3
t2 = t.*t;
t3 = t2.*t;
a(:,1) = 1.0;
a(:,2) = t;
a(:,3) = t2;
a(:,4) = t3;
case 4
t2 = t.*t;
t3 = t2.*t;
t4 = t3.*t;
a(:,1) = 1.0;
a(:,2) = t;
a(:,3) = t2;
a(:,4) = t3;
a(:,5) = t4;
case 5
t2 = t.*t;
t3 = t2.*t;
t4 = t3.*t;
t5 = t4.*t;
a(:,1) = 1.0;
a(:,2) = t;
a(:,3) = t2;
a(:,4) = t3;
a(:,5) = t4;
a(:,6) = t5;
case 6
t2 = t.*t;
t3 = t2.*t;
t4 = t3.*t;
t5 = t4.*t;
t6 = t5.*t;
a(:,1) = 1.0;
a(:,2) = t;
a(:,3) = t2;
a(:,4) = t3;
a(:,5) = t4;
a(:,6) = t5;
a(:,7) = t6;
case 7
t2 = t.*t;
t3 = t2.*t;
t4 = t3.*t;
t5 = t4.*t;
t6 = t5.*t;
t7 = t6.*t;
a(:,1) = 1.0;
a(:,2) = t;
a(:,3) = t2;
a(:,4) = t3;
a(:,5) = t4;
a(:,6) = t5;
a(:,7) = t6;
a(:,8) = t7;
otherwise
disp('Illegal polynomial degree')
end
col = degree+2;
for i = 1:nharmonics
iw = i*w;
a(:,col) = cos(iw.*t);
col = col+1;
a(:,col) = sin(iw.*t);
col = col+1;
end
y = a*coefs;
% Check derivatives flag
if ~derivs
fit = y;
return;
end
% Calculate yp
ap = zeros(npts,ncoefs);
switch degree
case 0
ap(:,1) = 0.0;
case 1
ap(:,1) = 0.0;
ap(:,2) = 1;
case 2
ap(:,1) = 0.0;
ap(:,2) = 1;
ap(:,3) = 2*t;
case 3
ap(:,1) = 0.0;
ap(:,2) = 1;
ap(:,3) = 2*t;
ap(:,4) = 3*t2;
case 4
ap(:,1) = 0.0;
ap(:,2) = 1;
ap(:,3) = 2*t;
ap(:,4) = 3*t2;
ap(:,5) = 4*t3;
case 5
ap(:,1) = 0.0;
ap(:,2) = 1;
ap(:,3) = 2*t;
ap(:,4) = 3*t2;
ap(:,5) = 4*t3;
ap(:,6) = 5*t4;
case 6
ap(:,1) = 0.0;
ap(:,2) = 1;
ap(:,3) = 2*t;
ap(:,4) = 3*t2;
ap(:,5) = 4*t3;
ap(:,6) = 5*t4;
ap(:,7) = 6*t5;
case 7
ap(:,1) = 0.0;
ap(:,2) = 1;
ap(:,3) = 2*t;
ap(:,4) = 3*t2;
ap(:,5) = 4*t3;
ap(:,6) = 5*t4;
ap(:,7) = 6*t5;
ap(:,8) = 7*t6;
otherwise
disp('Illegal polynomial degree')
end
col = degree+2;
for i = 1:nharmonics
iw = i*w;
ap(:,col) = -iw*sin(iw.*t);
col = col+1;
ap(:,col) = iw*cos(iw.*t);
col = col+1;
end
yp = ap*coefs;
% Calculate ypp
app = zeros(npts,ncoefs);
switch degree
case 0
app(:,1) = 0.0;
case 1
app(:,1) = 0.0;
app(:,2) = 0.0;
case 2
app(:,1) = 0.0;
app(:,2) = 0.0;
app(:,3) = 2;
case 3
app(:,1) = 0.0;
app(:,2) = 0.0;
app(:,3) = 2;
app(:,4) = 6*t;
case 4
app(:,1) = 0.0;
app(:,2) = 0.0;
app(:,3) = 2;
app(:,4) = 6*t;
app(:,5) = 12*t2;
case 5
app(:,1) = 0.0;
app(:,2) = 0.0;
app(:,3) = 2;
app(:,4) = 6*t;
app(:,5) = 12*t2;
app(:,6) = 20*t3;
case 6
app(:,1) = 0.0;
app(:,2) = 0.0;
app(:,3) = 2;
app(:,4) = 6*t;
app(:,5) = 12*t2;
app(:,6) = 20*t3;
app(:,7) = 30*t3;
case 7
app(:,1) = 0.0;
app(:,2) = 0.0;
app(:,3) = 2;
app(:,4) = 6*t;
app(:,5) = 12*t2;
app(:,6) = 20*t3;
app(:,7) = 30*t3;
app(:,8) = 42*t3;
otherwise
disp('Illegal polynomial degree')
end
col = degree+2;
w2 = w*w;
for i = 1:nharmonics
iw = i*w;
i2 = i*i;
i2w2 = i2*w2;
app(:,col) = -i2w2*cos(iw.*t);
col = col+1;
app(:,col) = -i2w2*sin(iw.*t);
col = col+1;
end
ypp = app*coefs;
% Store results
fit = [y yp ypp];
end