%% OpenAndFixFile
% Open and read the csv file. The Raven software will let you export the
% Raven data as a csv file. It's formatted incorrectly, but you can read
% it. Change the file below to your actual Raven .csv file. Most of this
% file is used to plot the power spectral density of the acceleromenter
% data but you can use the principles to read any of the Raven's sensors
fid  = fopen('E190_9SEP2017_J825R_1_Raven.csv','r');
f=fread(fid,'*char')';
fclose(fid);

% Fix the ,,-to-,, problem and write temporary file
goodFile = replace(f,",,",",,,");
fid  = fopen('temp.csv','w');
fprintf(fid,'%s',goodFile);
fclose(fid);

% Read the temporary file into a table
FIP_Data=readtable('temp.csv');
clearvars f goodFile fid

%% GetAccelData
timeAccel=transpose(FIP_Data.Time_AxialAccel_Gs_(isfinite(FIP_Data.Time_AxialAccel_Gs_)));
dtA = mean(diff(timeAccel))
accelZ=transpose(FIP_Data.AxialAccel_Gs_(isfinite(FIP_Data.AxialAccel_Gs_)));
timeAccelLat=transpose(FIP_Data.Time_LateralAccel_Gs_(isfinite(FIP_Data.Time_LateralAccel_Gs_)));
dtAL = mean(diff(timeAccelLat))
accelX=transpose(FIP_Data.LateralAccel_Gs_(isfinite(FIP_Data.LateralAccel_Gs_)));
plot(timeAccel,accelZ,timeAccelLat,accelX)
title('Accelerometer Data')
xlabel('time (s)')
ylabel('Acceleration')
legend('a_z', 'a_x')

%% Set Apogee Time
apogeeTime = 2;
timeAccelBoost = timeAccel(timeAccel<apogeeTime);
accelZBoost = accelZ(timeAccel<apogeeTime);
timeAccelBoostLat = timeAccelLat(timeAccelLat<apogeeTime);
accelXBoost = accelX(timeAccelBoostLat<apogeeTime);
plot(timeAccelBoost,accelZBoost,timeAccelBoostLat,accelXBoost)
title('Accelerometer Data during Boost')
xlabel('time (s)')
ylabel('Acceleration')
legend('a_z', 'a_x')

%% RemoveBaseline
accelZbl = accelZ - mean(accelZ);
accelXbl = accelX - mean(accelX);
plot(timeAccel,accelZbl,timeAccelLat,accelXbl)
title('Accelerometer Data')
xlabel('time (s)')
ylabel('Acceleration')
legend('a_z', 'a_x')

%% DoPSD
[PSDX, fx] = pwelch(accelXbl,[],[],[],1/dtAL);
[PSDZ, fz] = pwelch(accelZbl,[],[],[],1/dtA);
semilogx(fx,10*log10(PSDX),fz,10*log10(PSDZ))
title('Welch PSD Plots')
xlabel('Frequency (Hz)')
ylabel('Acceleration in dB (g/Hz)')
legend('a_x', 'a_z')

%% RemoveBaseline
accelZblB = accelZBoost - mean(accelZBoost);
accelXblB = accelXBoost - mean(accelXBoost);
plot(timeAccelBoost,accelZblB,timeAccelBoostLat,accelXblB)
title('Accelerometer Data')
xlabel('time (s)')
ylabel('Acceleration')
legend('a_z', 'a_x')

%% DoPSD
[PSDX, fx] = pwelch(accelXblB,[],[],[],1/dtAL);
[PSDZ, fz] = pwelch(accelZblB,[],[],[],1/dtA);
semilogx(fx,10*log10(PSDX),fz,10*log10(PSDZ))
title('Welch PSD Plots')
xlabel('Frequency (Hz)')
ylabel('Acceleration in dB (g/Hz)')
legend('a_x', 'a_z')


