% This script tests what happens when one of the sinusoidal components of 
% the system does not have data for an entire period.

% Number of time series
num_series = 6;

% Interval end points
t_low = 0;
t_high = 2.5;
t_high = 16;

% Time step size
del_t = 0.1;
%del_t = 0.01;

% Time Points
T = [t_low:del_t:t_high];

% Number of time points
N = length(T);

% Amplitudes, frequencies, and phase shifts of components
a1 = 1;
w1 = 2*pi;
phi1 = 2*pi*rand(num_series,1);
a2 = 10;
w2 = 0.25*pi;
phi2 = 2*pi*rand(num_series,1);

% Construct signal and remove the mean for each data series
X = a1*sin(ones(num_series,1)*w1*T+phi1*ones(1,N)) +  ...
    a2*sin(ones(num_series,1)*w2*T+phi2*ones(1,N));
X = X - (mean(X'))'*ones(1,N);

% SVD of data set
[u,s,v] = svd(X);
s = diag(s)

% Plot time series for first 4 variables
figure(1);
hold off;
plot(T,X(1,:));
hold on;
plot(T,X(2,:));
plot(T,X(3,:));
plot(T,X(4,:));

hold off;

% Plot first 4 "modes"
figure(2);
plot(T',v(:,1));
hold on;
plot(T',v(:,2));
plot(T',v(:,3));
plot(T',v(:,4));
%plot(T',v(:,5));
hold off;

% Do a fourier analysis of the first few "modes"
max_freq = 20;
W = [0:max_freq-1]*2*pi/N/del_t;
F1 = fft(v(:,1));
F2 = fft(v(:,2));
F3 = fft(v(:,3));
F4 = fft(v(:,4));
figure(3);
plot(W,abs(F1(1:max_freq)));
hold on;
plot(W,abs(F2(1:max_freq)));
plot(W,abs(F3(1:max_freq)));
plot(W,abs(F4(1:max_freq)));
hold off;