%% Typical Shock Accel Processing
% This example demonstrates the typical steps for analyzing accelerometer
% data physically measured in a drop test of a consumer electronic device. 
% The example takes the user through the key steps of data analysis:
% assessing the raw data, performing several data quality checks using
% Fourier spectrum analysis and integrating acceleration to velocity, 
% applying lowpass filtering to improve data interpretation, and finally
% computing a SRS for the data. 

%% Set-up
%{
This section makes a few initial settings for the example.
 
Additional Details: 
The settings below include: 
 o Popup asking to clean MATLAB session before start of example.
 o Define a handy "adv" variable so that we can later use it to 
   quickly access any Kornucopia function's ADV options with ease.
 o Setting a specific Units Preference for this example. 
 o Define some units variables that will be used throughout the 
   example.
 o Turn off a warning regarding Props not being inheritted during
   certain math operations.
 o Create a tabbed figure with the main tab titles that will be used for
   the entire analysis
%}

k_cleanup();
adv = k_adv();
k_unitsPreferenceActivate('ShockVib_mm_ms_G');
k_unitsVariables('ms, G, Hz, kHz');
k_set('inheritPropsAutoWarning', 'off');

figH = figure('Name', 'Typical Shock Accel Processing', ...
    'pos', k_figZoom(2));

tabStuffMain = k_figTabsCreate(6, figH);

txt = 'Run the next section in the script to start filling in the figure.';
uiwait(msgbox(txt, 'Notice'), 10);

%% Raw Data Assessment

% Load raw accel data
dataDir = '..\data';
fileName = 'accel_data.txt'; % Use string '?' for file selection popup
fName = fullfile(dataDir, fileName);
delim = ',';
startStr = [];
rs1 = [];
colHeadStr = 'Time';
rs2 = 0;
endStr = [];
rs3 = [];
raw = k_readText(fName,delim,startStr,rs1,colHeadStr,rs2,endStr,rs3);

% Retrieve descriptive header info from the file.
header = raw.Props.UserData.MainFileHeader;

% Replace the variable's description with the Header info 
raw.Props.Description = header;

% Compute sampling Freq of data
fs = k_samplingFreq(raw('Time'));
fs = fs.convert('MHz');

% Present the raw data
currentTabH = tabStuffMain.TabsH(1);
k_figTabsClear(currentTabH);
currentTabH.Title = 'Raw Accel Plots';

k_displayOnFigure(0.25, 0.75, raw, ...
    'parent', currentTabH, ...
    'box', 'on', ...
    'removeVarName', 'off')
    
tmpH = cell(3,1);
toPlt = raw;

subplot(2,2,2, 'parent', currentTabH);
tmpH{1} = k_plot(toPlt, [], ...
    'depend', 'X', ...
    'grid', 'major');

txt = ['Peak abs Accel = ', num2str(max(abs(toPlt('X'))),3)];
k_displayOnFigure(0.85, 0.9, txt, ...
    'box', 'on');
    
subplot(2,2,3);
tmpH{2} = k_plot(toPlt, [], ...
    'depend', 'Y', ...
    'grid', 'major');

txt = ['Peak abs Accel = ', num2str(max(abs(toPlt('Y'))),3)];
k_displayOnFigure(0.85, 0.9, txt, ...
    'box', 'on');

subplot(2,2,4);
tmpH{3} = k_plot(toPlt, [], ...
    'depend', 'Z', ...
    'grid', 'major');

txt = ['Peak abs Accel = ', num2str(max(abs(toPlt('Z'))),3)];
k_displayOnFigure(0.85, 0.9, txt, ...
    'box', 'on');

k_displayOnFigure(0.5, 0.97, fs, ...
    'parent', currentTabH, ...
    'removeVarName', 'off', ...
    'fontSize', 16, ...
    'fontWeight', 'Bold');

k_plotMatchLimits(tmpH, 'xy');

k_figTabsDisplay(currentTabH);

%% Using Pre-trigger Section, Compute and Remove DC Amplitude Offsets

%{
The function k_offsetAdjustY removes undesirable amplitude offsets in 
data by averaging (by column) the yData in a pre-trigger region (specified
by the "offsetRegionEndVal" argument) and then subtracting these average 
offsets from all the data in a given column.

We make this adjust so that during integration to assess velocity we do
not get accessive distortion due to an initial DC bias (offset) in the
accel data.
%}

adj = k_offsetAdjustY(raw, [], '-1*ms');

% Plot the adjusted data
currentTabH = tabStuffMain.TabsH(2);
k_figTabsClear(currentTabH);
currentTabH.Title = 'After Removing DC Bias';

txt = {
    'Data adjusted to remove mean DC amplitude offsets.'
    'Small changes in Peak Accels occur from this adjustment.'
    ' '
    'Note: The data at this point is still essentially "raw", '
    'meaning it contains the intended signals PLUS potentially '
    'some distortion from the accelerometer''s transfer function. ' 
    };
k_displayOnFigure(0.25, 0.75, txt, ...
    'parent', currentTabH, ...
    'box', 'on', ...
    'fontSize', 12)

toPlt = adj;
cols = {'X', 'Y', 'Z'};
for q = 1 : 3
    subplot(2,2,q+1, 'parent', currentTabH);
    tmpH{q} = k_plot(toPlt, [], ...
        'depend', cols{q}, ...
        'grid', 'major');
    
    txt = ['Peak abs Accel = ', num2str(max(abs(toPlt(cols{q}))),3)];
    k_displayOnFigure(0.85, 0.9, txt, ...
        'box', 'on');    
end   

k_plotMatchLimits(tmpH, 'xy');

k_figTabsDisplay(currentTabH);

%% Compute DFS of Data

dfs = k_DFS(adj, [], 'hann', 'mag');

% Plot the results
currentTabH = tabStuffMain.TabsH(3);
k_figTabsClear(currentTabH);
currentTabH.Title = 'DFS of Data';

txt = {
    'From the DFS plots we can clearly see the '
    'resonance of the Endevco 727 accelerometers. '
    'This is common for MEMS-based piezoresistive sensors.'
    ' '
    'For the 3 accels, their resonance is approximately '
    'fnSensor = 200*kHz or higher. '
    ' '
    'The max bandwidth of the measurement that would be'
    'undistorted by the sensor''s transfer function'
    'resonance behavior is 0.1*fnSensor = 20*kHz.'    
    'This is presented on the graphs by the Orange line.'
    };
k_displayOnFigure(0.25, 0.75, txt, ...
    'parent', currentTabH, ...
    'box', 'on', ...
    'fontSize', 12)

fnSensor = 200*kHz;
fc = 0.1*fnSensor;

toPlt = dfs;
for q = 1 : 3
    subplot(2,2,q+1, 'parent', currentTabH);
    tmpH{q} = k_plot(toPlt, [], ...
        'depend', cols{q}, ...
        'scale', 'logLog', ...
        'xGuideVal', fc);
end   

k_plotMatchLimits(tmpH, 'xy');

k_figTabsDisplay(currentTabH);

%% Lowpass Filter the Data

% Compute filter coefficents
c = k_cIIR(fs, fc, 'lp', 'butter'); % default filter order is 6
filtered = k_filter(adj, [], c, 2); % last arg, 2, implies bi-directional

% Plot the results
currentTabH = tabStuffMain.TabsH(4);
k_figTabsClear(currentTabH);
currentTabH.Title = 'Lowpass Filtering to 0.1*fnSensor';

txt = {
    'These are the best representations of the accelerations'
    'that occurred in the product being tested. '
    ' '
    'Notes: '
    '1) The Y and Z sensors now show much lower response. '
    '   What we were seeing before was mostly sensor resonance.'
    '2) Look back at the Raw Accel tab (or Adjusted data tab).'
    '   Those results were being distorted by the sensor''s '
    '   transfer function resonance. '
    };
k_displayOnFigure(0.05, 0.75, txt, ...
    'parent', currentTabH, ...
    'box', 'on', ...
    'alignHor', 'left', ...
    'fontSize', 12);

toPlt = filtered;
for q = 1 : 3
    subplot(2,2,q+1, 'parent', currentTabH);
    tmpH{q} = k_plot(toPlt, [], ...
        'depend', cols{q}, ...
        'grid', 'major');
    
    txt = ['Peak abs Accel = ', num2str(max(abs(toPlt(cols{q}))),3)];
    k_displayOnFigure(0.85, 0.9, txt, ...
        'box', 'on');    
end   

k_plotMatchLimits(tmpH, 'xy');

k_figTabsDisplay(currentTabH);

%% Velocity Checks

%{
We integrate the acceleration to obtain velocity vs time. We can use this 
as a sanity check of the acceleration data. 

For free-drop testing, the velocity change we measure is expected to be
bounded by the relationship between potential energy for a drop of height 
h, PE=(M*G*h) and the kinetic energy just prior to impact KE=(0.5*M*v^2). 

Initial velocity prior to impact = sqrt(2*G*h);
Thus, the deltaV in the direction of impact should be bounded by this 
initial impact velocity (perfectly plastic impact) and twice this velocity
(perfectly elastic impact).

%}

% Extract the drop height from the last row of the variable's description
hDrop = k_unpackHeader(raw.Props.Description{end}, '=');

% Compute deltaV theoretical bounds
deltaV_theoryBounds = [sqrt(2*G*hDrop), 2*sqrt(2*G*hDrop)];

% Integrate accel to obtain transient velocities
Vel = k_integrate(filtered);

% Plot the results
currentTabH = tabStuffMain.TabsH(5);
k_figTabsClear(currentTabH);
currentTabH.Title = 'Velocity Checks';

toPlt = Vel;
k_plot(toPlt, [], ...,
    'parent', currentTabH, ...
    'grid', 'major', ...
    'yGuideVals', deltaV_theoryBounds, ...
    'yLim', '-inf*mm/ms, 8*mm/ms', ...
    'yLabelBase', {'user', 'Velocity'})

k_displayOnFigure(0.2, deltaV_theoryBounds(1,1), 'Lower Bound, CR = 0' , ...
    adv.k_displayOnFigure.locationUnits.mixedNormAndData, ...
    'fontWeight', 'bold', ...
    'box', 'on', ...
    'boxLineStyle', 'none');

k_displayOnFigure(0.2, deltaV_theoryBounds(1,2), 'Upper Bound, CR = 1' , ...
    adv.k_displayOnFigure.locationUnits.mixedNormAndData, ...
    'fontWeight', 'bold', ...
    'box', 'on', ...
    'boxLineStyle', 'none');


txt = {
    'Velocity change of X-direction sensor falls within '
    'theoretical bounds for coefficient of restitution (CR).'
    };
k_displayOnFigure(0.75, 0.75, txt, ...
    'box', 'on', ...
    'fontSize', 12);

txt = {
    'Sensors Y and Z are not in primary drop direction,'
    'so their velocity changes are expected to be relatively small.'
    };
k_displayOnFigure(0.75, 0.28, txt, ...
    'box', 'on', ...
    'fontSize', 12);

k_figTabsDisplay(currentTabH);

%% Shock Response Spectrum

%{
In this section we compute the SRS on both the adjusted data
and the filtered data. 
   o You can specify the SRS type as 'PVSS' or 'Accel'
   o Another common SRS type of calculation would be PVSS. If you choose 
     to compute a PVSS you should consider using the following ADV options
     for k_plot:
        'scale', 'fourCP', ... 
        'yLabelBase', {'user', 'Pseudo Vel'}        
%}

% Some SRS parameters
fnMin = 10*Hz;
fnMax = 'fnMax';
damping = 0.05;
SRStype = 'PVSS'; % Options are 'Accel' or 'PVSS'

srs = struct;

% SRS on adjusted data
srs.adj = k_SRS(adj, [], fnMin, fnMax, damping, SRStype, ...
    'modifyColNames', 'replace');

srs.adj.Props.Description = {
    'SRS on Data Including Sensor Resonance'
    srs.adj.Props.Description{end}
    };

% SRS on LP filtered data
srs.filtered = k_SRS(filtered, [], fnMin, fnMax, damping, SRStype, ...
    'modifyColNames', 'replace');

srs.filtered.Props.Description = {
    'SRS on Filtered Data (fc = 0.1*fnSensor)'
    srs.filtered.Props.Description{end}
    };

% Plot the results
currentTabH = tabStuffMain.TabsH(6);
k_figTabsClear(currentTabH);
currentTabH.Title = 'Shock Response Spectrum';

tmpH = cell(2,1);

if strcmpi(SRStype, 'Accel')
    SRSpltADVs = {'scale', 'logLog', ...
        'xGuideVal', fc, ...
        'lineWidth', 2, ...
        'title', 'description'};    
else
    % PVSS
    SRSpltADVs = {'scale', 'fourCP', ...
        'xGuideVal', fc, ...
        'lineWidth', 2, ...        
        'title', 'description', ...
        'legendLocation', 'SouthWest'};     
end

subplot(2,2,1, 'parent', currentTabH);
tmpH{1} = k_plot(srs.adj, [], ...
    SRSpltADVs);

subplot(2,2,2);
tmpH{2} = k_plot(srs.filtered, [], ...
    SRSpltADVs);

txt = {
    'The orange line in plots above denotes fc.'
    ' '
    'By filtering the data with fc = 0.1*fnSensor for the'
    'plot on the right, we removed distoritions caused by the '
    'sensor resonance and achieve the best estimate of the SRS. '
    };
k_displayOnFigure(0.5, 0.45, txt, ...
    'parent', currentTabH, ...
    'box', 'on', ...
    'alignVert', 'top', ...
    'fontSize', 12);

k_plotMatchLimits(tmpH, 'xy');

k_figTabsDisplay(currentTabH);

%% Options to save figure windows

k_exampleSetup('figSavingOptions');
