% Example analysis of vibration data.
% The accompanying s2_tailoff.txt is accelerometer data from an ACS anomaly during
% stage 2 tailoff for a rocket. The data are from Tom Irvine, and are used by 
% permission. For many more good examples, see his website, vibrationdata.com. 


% Get the data and plot it.
M = dlmread("s2_tailoff.txt");
t = M(:,1);          % time in seconds
dt = t(2) - t(1);
accel = M(:,2);      % acceleration in g

figure(1);
hold off
plot(t, accel);
hold on
xlabel('t (s)')
ylabel('acceleration (g)')
title('Stage 2 Tailoff: Accelerometer Data')
print('../plots/s2_tailoff_accel.pdf','-dpdf')
axis([124,125])
print('../plots/s2_tailoff_accel_1s.pdf','-dpdf')
axis([125,125.1])
print('../plots/s2_tailoff_accel_0.1s.pdf','-dpdf')


% Region of interest.
roi = find(t>124.71 & t<131.55); % Define region of interest.
a_roi = accel(roi);

% Power spectrum (spectrum.m requires a binary file of big-endian floats).
Nchunk = 256;  % chunk size
file_handle = fopen('s2_tailoff.bin','w');
fwrite(file_handle, a_roi, "float", "ieee-be");
fclose('s2_tailoff.bin'); % *** This is rather important ***
[powerspec,frequency] = spectrum('s2_tailoff.bin', Nchunk, dt);
df = frequency(2) - frequency(1);
rms_accel = std(accel);  % calculate RMS acceleration from data
normalization = rms_accel / sum(powerspec)*df;  
   % take advantage of Parseval's theorem no normalize the spectrum.
powerspec = powerspec * normalization;  % units: g^2 / Hz

figure(2);
hold off
loglog(frequency(2:end),powerspec(2:end))
hold on
title('Stage 2 Tailoff: PSD')
xlabel('frequency (Hz)')
ylabel('power (g^2/Hz)')
print('../plots/s2_tailoff_PSD.pdf','-dpdf')



% Vibration kinematics
% There is an important period at around 40 ms (25 Hz), which is
% readily apparent in the envelope of the accelerometer data.
% Optimize the filter to explore the structure within that envelope: 
Nfilt = 0.04/dt; 

velocity = ifilter(accel' .* 9.81 .* dt, Nfilt);
   % Notes: 
   % (1) converting accel to m/s^2.
   % (2) took awhile to realize that I had to transpose accel...!
velocity_unfilt = cumsum(accel')*9.81*dt;

figure(3);
hold off
plot(t, velocity, 'b')
hold on
title('Stage 2 Tailoff: Velocity')
xlabel('t (s)')
ylabel('velocity (m/s)')
legend('velocity (40 ms filter)')
print('../plots/s2_tailoff_vel.pdf','-dpdf')
plot(t, velocity_unfilt, 'r')
legend('velocity (40 ms filter)', 'velocity (unfiltered)')
print('../plots/s2_tailoff_vel_unfilt.pdf','-dpdf')
axis([124, 125, -0.3, 0.1])
print('../plots/s2_tailoff_vel_1s.pdf','-dpdf')
axis([125, 125.1, -0.05, 0.05])
legend('velocity (40 ms filter)')  % unfiltered not visible at this zoom.
print('../plots/s2_tailoff_vel_0.1s.pdf','-dpdf')

displacement = ifilter(velocity*dt, Nfilt);
displacement_unfilt = cumsum(velocity_unfilt)*dt;

figure(4);
hold off
plot(t, displacement, 'b')
hold on
title('Stage 2 Tailoff: displacement')
xlabel('t (s)')
ylabel('displacement (m)')
legend('displacement (40 ms filter)')
print('../plots/s2_tailoff_disp.pdf','-dpdf')
plot(t, displacement_unfilt, 'r')
legend('displacement (40 ms filter)','displacement (unfiltered)')
print('../plots/s2_tailoff_disp_unfilt.pdf','-dpdf')
axis([124, 125, -0.003, 0.001])
print('../plots/s2_tailoff_disp_1s.pdf','-dpdf')
axis([125, 125.1, -0.0001, 0.0001])
legend('displacement (40 ms filter)')  % unfiltered not visible at this zoom.
print('../plots/s2_tailoff_disp_0.1s.pdf','-dpdf')


% Some printed results --- statistics on vibrations within the ROI:
roi_time =  t(roi(end)) - t(roi(1))
rms_accel = std(accel(roi))
rms_velocity = std(velocity(roi))
rms_displacement = std(displacement(roi))
total_travel = sum(abs(displacement(roi)))
mean_speed = total_travel/roi_time