function [ ] = nickFunction01( dataCell, trimPointsStart, movingAverageSize )

% step 0: input cellArray of IMU data with all rows using the Matlab inport tool
% step 1: convert to matrix of floats
% step 2: detrend the data to remove DC offset
% step 2: obtain welch spectrum, plot both time and freq domains
% step 3: apply low pass filter, plot both time and freq domains
% step 4: integrate the filtered data to obtain velocity estimate
close all

% Input data format
% 1. sample #
% 2. accel x 
% 3. accel y
% 4. accel z
% 5. ...

%   Detailed explanation goes here

%% Step 1: convert cellArray to floatArray

% the number of columns to read from cellarray
readColumns = [3,4,5,19,20,21];

% convert from cell array to matrix
dataImport = cell2mat(dataCell(trimPointsStart:end, readColumns));
numSamples = linspace(0,size(dataImport,1),size(dataImport,1))';

% [sample#, accel x, accel y, accel z....]
dataFloat = [numSamples, dataImport];

accel = zeros(size(dataImport,1),3);

% convert from m/sec^2 into mm/sec^2
for i = 1:3
    accel(:,i) = dataFloat(:,2+i).*-1000;
end

figure(1)
plot(numSamples,accel(:,1), 'r-', numSamples, accel(:,2), 'g-', numSamples, accel(:,3),'b-')
title('Raw Sensor Data')
xlabel('Samples');
ylabel('Acceleration Magnitude (mm/s^2)');
legend('X accel', 'Y accel', 'Z accel', 'best');


%% Step 2: detrend the data
accel_detrend = zeros(size(dataImport,1),3);

for i = 1:3
    accel_detrend(:,i) = detrend(accel(:,i));
end

figure(2)
plot(numSamples,accel_detrend(:,1), 'r-', numSamples, accel_detrend(:,2), 'g-', numSamples, accel_detrend(:,3),'b-')
title('Raw Sensor Data')
xlabel('Samples');
ylabel('Acceleration Magnitude (mm/s^2)');
legend('X accel', 'Y accel', 'Z accel', 'best');


%% Step 2: obtain welch spectrum

% sampling freq of signal
fs = 100; %[Hz]

% length of DFT
nfft = 256;

% time values form input signal
time = 1/fs.*numSamples; % [sec]

% length of window 
window = hamming(nfft);

% number of overlapped samples
noverlap = length(window)/2;

% PDS of raw data
[accel_welch,f] = pwelch(accel, window, noverlap, nfft, fs);

% PDS of detrended data
[accel_detrend_welch,f_detrend] = pwelch(accel_detrend, window, noverlap, nfft, fs);

figure(3)
plot(f, 10*log10(accel_welch(:,1)), 'r-', f, 10*log10(accel_welch(:,2)), 'g-', f, 10*log10(accel_welch(:,3)), 'b-')
title('Welch Method - Raw Data');
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
legend('X axis', 'Y axis', 'Z axis', 'best');
grid

figure(4)
plot(f_detrend, 10*log10(accel_detrend_welch(:,1)), 'r-', f_detrend, 10*log10(accel_detrend_welch(:,2)), 'g-', f_detrend, 10*log10(accel_detrend_welch(:,3)), 'b-')
title('Welch Method - Raw Data');
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
legend('X axis', 'Y axis', 'Z axis', 'best');
grid


%% Step 03: Filter data

% moving average filter
a = 1;
b = ones(1,movingAverageSize)./movingAverageSize;

% Frequency response of filter
[H,F] = freqz(b,a,1024,fs);

figure(10)
plot(F,10*log10(H),'-');
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
title('Frequency Response of Moving Average Filter');

% Apply moving average filter to detrend accel data
accel_detrend_filtered = filter(b,a,accel_detrend);
accel_detrend_filtfilt = filtfilt(b,a,accel_detrend);


figure(5)
hold on
plot(numSamples,accel_detrend_filtered(:,1), 'r-', numSamples, accel_detrend_filtered(:,2), 'g-', numSamples, accel_detrend_filtered(:,3),'b-')
%plot(numSamples,accel_detrend_filtfilt(:,1), 'r--', numSamples, accel_detrend_filtfilt(:,2), 'g--', numSamples, accel_detrend_filtfilt(:,3),'b--')
hold off
title('Filtered Sensor Data')
xlabel('Samples');
ylabel('Acceleration Magnitude (mm/s^2)');
legend('X accel', 'Y accel', 'Z accel', 'best');

%% Step 04: Numerical integration

%reference: http://www.mathworks.com/help/matlab/ref/cumtrapz.html

velocity = cumtrapz(accel_detrend).*(1/fs);
velocity_filtered = cumtrapz(accel_detrend_filtered).*(1/fs);

figure(6)
plot(time, velocity(:,1), 'r-',time, velocity(:,2), 'g-', time, velocity(:,3), 'b-')
title('Numerical Integration - Trapezoid Rule');
xlabel('Time (sec)')
ylabel('Velocity (mm/s)')
legend('X velocity', 'Y velocity', 'Z velocity', 'best');

figure(7)
plot(time, velocity_filtered(:,1), 'r-',time, velocity_filtered(:,2), 'g-', time, velocity_filtered(:,3), 'b-')
title('Numerical Integration Filtered - Trapezoid Rule');
xlabel('Time (sec)')
ylabel('Velocity (mm/s)')
legend('X velocity', 'Y velocity', 'Z velocity', 'best');



%% Step 05: total magnitude 

% calculate the magnitude from individual velocity estimates
velocity_magnitude = sqrt(velocity(:,1).^2 + velocity(:,2).^2 + velocity(:,3).^2);
velocity_magnitude_filtered = sqrt(velocity_filtered(:,1).^2 + velocity_filtered(:,2).^2 + velocity_filtered(:,3).^2);

% calculate the magnitude from filtered acceleration data
accel_magnitude = sqrt(accel_detrend_filtered(:,1).^2 + accel_detrend_filtered(:,2).^2 + accel_detrend_filtered(:,3).^2);
velocity_estimate = cumtrapz(accel_magnitude).*(1/fs);

% calculate the magnitude from acceleration detrended data
accel_detrend_magnitude = sqrt(accel_detrend_filtered(:,1).^2 + accel_detrend_filtered(:,2).^2 + accel_detrend_filtered(:,3).^2);
velocity_estimate = cumtrapz(accel_detrend_magnitude).*(1/fs);
velocity_estimate_detrend = detrend(velocity_estimate);

figure(8)
%plot(time, velocity_magnitude, 'r-', time, velocity_magnitude_filtered, 'g-', time, velocity_estimate, 'b-', time, velocity_estimate_detrend)
plot(time, velocity_magnitude, 'r-', time, velocity_magnitude_filtered, 'b-')

title('Velocity Magnitude');
xlabel('Time (sec)')
ylabel('Velocity (mm/s)')
legend('Velocity - raw', 'Velocity - filtered','Velocity - integrate from accel', 'Velocity - integrate from accel detrend', 'best');


end