function [ ] = nickFunction02( dataCell, trimPoints, movingAverageSize )

% Description: Input both the raw accelerometer and IMU estiamte acceleration data

% 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

% number of points to remove trim from start of data set
%trimPoints = 150;

% conver from cell array to matrix
dataImport = cell2mat(dataCell(trimPoints:end, 1:6));
numSamples = linspace(0,size(dataImport,1),size(dataImport,1))';

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

% convert from m/sec^2 into mm/sec^2
accel_raw = dataFloat(:,3).*-1000;
accel_IMU = dataFloat(:,6).*-1000;

figure(1)
hold on
plot(accel_raw,'b-')
plot(accel_IMU, 'g--')
hold off
title('Raw Data')
xlabel('Samples');
ylabel('Acceleration Magnitude (mm/s^2)');
legend('raw data', 'IMY data', 'Location', 'SouthWest');

%% Step 2: detrend the data
accel_raw_detrend = detrend(accel_raw);
accel_IMU_detrend = detrend(accel_IMU);

figure(2)
hold on
plot(accel_raw_detrend,'b-')
plot(accel_IMU_detrend, 'g--');
hold off
title('After Detrending');
xlabel('Samples');
ylabel('Acceleration Magnitude (mm/s^2)');
legend('raw data', 'IMU data', 'Location', 'SouthWest');



%% 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;

% Determine spectrum power density (no detrend)
[pxx_raw,f_raw] = pwelch(accel_raw, window, noverlap, nfft, fs);
[pxx_IMU,f_IMU] = pwelch(accel_raw, window, noverlap, nfft, fs);

% Determine spectrum power density (with detrending)
[pxx_raw_detrend,f_raw_detrend] = pwelch(accel_raw_detrend, window, noverlap, nfft, fs);
[pxx_IMU_detrend,f_IMU_detrend] = pwelch(accel_raw_detrend, window, noverlap, nfft, fs);

figure(3)
hold on
plot(f_raw, 10*log10(pxx_raw), 'b-')
plot(f_IMU, 10*log10(pxx_IMU), 'g--')
hold off
title('Welch Method - (no detrend)');
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
legend('raw data (no detrend)', 'IMU data (no detrend)', 'Location', 'SouthWest');
grid

figure(4)
hold on
plot(f_raw_detrend, 10*log10(pxx_raw_detrend),'b-')
plot(f_IMU_detrend, 10*log10(pxx_IMU_detrend),'g--')
hold off
title('Welch Method - (with detrend)');
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
legend('raw data (with detrend)', 'IMY data (with detrend)', 'Location', 'SouthWest');
grid

%% Step 03: Filter data

% moving average filter
a = 1;
b = ones(1,movingAverageSize)./movingAverageSize;
[H,F] = freqz(b,a,1024,fs);

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

% Apply filter to data after detrending
y_raw_filter = filter(b,a,accel_raw_detrend);
y_IMU_filter = filter(b,a,accel_IMU_detrend);

figure(5)
hold on
plot(accel_raw_detrend, '-k');
plot(y_raw_filter, '--b');
hold off
xlabel('Samples');
ylabel('Acceleration Magnitude (mm/s^2)');
grid on
legend('Raw - Original Data','Filtered Data','Location','SouthWest')
title('Plot of Original and Filtered Data (Raw sensor output)')

figure(6)
hold on
plot(accel_IMU_detrend, '-k');
plot(y_IMU_filter, '--b');
hold off
xlabel('Samples');
ylabel('Acceleration Magnitude (mm/s^2)');
grid on
legend('IMU - Original Data','Filtered Data','Location','SouthWest')
title('Plot of Original and Filtered Data (IMU output')



%% Step 04: Numerical integration

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

Z_raw = cumtrapz(accel_raw_detrend).*(1/fs);
Y_raw = cumtrapz(y_raw_filter).*(1/fs);

Z_IMU = cumtrapz(accel_IMU_detrend).*(1/fs);
Y_IMU = cumtrapz(y_IMU_filter).*(1/fs);

figure(7)
hold on 
plot(time, Z_raw, 'k-');
plot(time, Y_raw, 'b--');
plot(time, Z_IMU, 'g-');
plot(time, Y_IMU, 'r--');
hold off
title('Numerical Integration - Trapezoid Rule');
xlabel('Time (sec)')
ylabel('Velocity (mm/s)')
legend('Raw - Original Data','Raw - Filtered Data','IMU - Original Data','IMU - Filtered Data','Location','SouthWest')


% figure(8)
% hold on 
% plot(time, Z_IMU, 'k-');
% plot(time, Y_IMU, 'b--');
% hold off
% title('IMU sensor output - Numerical Integration - Trapezoid Rule');
% xlabel('Time (sec)')
% ylabel('Velocity (mm/s)')
% legend('IMU - Original Data','IMU - Filtered Data','Location','SouthEast')



end