%% Prepare Workspace
clear; %clears variables
close all; %closes all figures

%% Prepare data
[file,path] = uigetfile('.csv'); %choose file
data = readtable(file); %get data from file in a table
time = data{:,1}; %extracting only the first column for time data
y = data{:,2}; %extracting only second column for the y data

%% FFT
nPoints = length(y); %number of samples
fs = 1500; %sample frequency

fullSpec = fft(y);
fullSpecpos = abs(fullSpec);
halfSpec = fullSpecpos(1:nPoints/2+1); 
halfSpec = halfSpec*2/nPoints; %normalise amplitude
halfSpec(1) = halfSpec(1)/2; %correct DC component

f = (fs*(0:nPoints/2))/nPoints; %normalise frequency

%% Interpolation
n = 150; %points for polynom;
degree = 2; %degree of polynom;
[maxA,idMax] = max(halfSpec);

%get three data points around the max and fit coefficients for poly
fPoly = [f(idMax-1) f(idMax) f(idMax+1)];
aPoly = [halfSpec(idMax-1) halfSpec(idMax) halfSpec(idMax+1)];
coeffs = polyfit(fPoly,aPoly,degree);
minfPoly = f(idMax-1);
maxfPoly = f(idMax+1);

%form poly
xPoly = linspace(minfPoly,maxfPoly,n);
yPoly = polyval(coeffs,xPoly);

%evaluate poly
dcoeffs = polyder(coeffs);
rootsPoly = roots(dcoeffs);
maxPolyA = polyval(coeffs,rootsPoly);


%% plot everything
%plot data
fig = figure();
subplot(2,1,1);
plot(time,y,'black');
xlim([0 0.1]);
title('Sampling Frequency fs = 1500Hz');
ylabel('y(t) [-]');
xlabel('Time t [s]');


%plot fft
subplot(2,1,2);
hold on; %to plot to functions in one plot
stem(f,halfSpec);
xlim([30 55]);
ylim([0 0.4]);
ylabel('Magnitude');
xlabel('Frequency f [Hz]');

%ploting the interpolation
plot(xPoly,yPoly,'red');
plot(minfPoly,halfSpec(idMax-1),'ro');
plot(f(idMax),halfSpec(idMax),'ro');
plot(f(idMax+1),halfSpec(idMax+1),'ro');
plot(rootsPoly,maxPolyA,'b.');
text(maxfPoly+2,maxPolyA,'fmax = '+string(rootsPoly) + ' Hz');
text(maxfPoly+2,maxPolyA-0.1,'a(fmax) = '+string(maxPolyA));
