%% Load Data

% The first column of the csv is time; the second is x(t).
data = readmatrix('test_data.csv');
t = data(:,1); x = data(:,2);

% Calculate "delta t", the sampling time. Obviously, this only works if
% the same delta t was used for each sample taken.
dt = t(2) - t(1);

%% Carry out FFT

% The total number of frequencies we can possibly investigate is equal to 
% N/2 + 1, where N is the length of the data vector avaialble to us.
N = length(x);

% Carry out the discrete Fourier Transform of this data, and find the
% magnitude of the resulting complex number.
X = fft(x);
mag = abs(X(1:N/2+1)) * (2/N);

% Maximum sampling rate is the inverse of the time step
sampling_rate = 1/dt;

% Create a vector of equally spaced "Hertz" values from 0 to sampling_rate
freqs =  2*pi * (0:N/2)' * (sampling_rate/N);

%% Make Plots

figure(2); clf;

subplot(3,1,1);
plot(t,x,'.-'); title("Data in time domain");
xlabel("Time [s]");

subplot(3,1,2);
plot(freqs,mag); title("Data in frequency domain");
xlabel("Angular Frequency [rad/s]")

%% Overlay identified function
% A new time axis has been created with 10 times the sampling frequency,
% with the goal of finding a function f(x) that reasonably captures the
% data using a small number of Fourier modes. The code below currently uses
% sin(x), which clearly does not capture the function well. 

% Identify the function by trial and error
figure(2); subplot(3,1,3);
plot(t,x,'.-'); title("Approximate Function"); hold on;
t_fine = linspace(0,t(end),10*length(t));
plot(t_fine,sin(t_fine)); % CHANGE THIS LINE ONLY
xlim([35,38]); xlabel("Time [s]")
legend("Data","Approximate Function");