%% Moon Penetrator
% This example demonstrates how to use DSP tools to correlate 
% physical test data and simulation data. The example presented here
% comes from evaluation of a penetrator intended to be sent to an 
% ice-covered moon. The primary data being assessed is acceleration vs 
% time, where the physical test and FEA simulations look nothing like 
% each other at the beginning of the analysis. By using various DSP tools, 
% in the end it is clear which simulation best correlates to the physical 
% test data. 

%% Overview
%{
Key concepts and features demonstrated include: 
 o Using a combination of DSP features within Kornucopia to assess and
   correlate the physical test data and the simulation results. 
   Features demonstrated include:
   - Efficiently reading in raw text files,
   - Regularizing the FEA data to a constant sampling increment,
   - Performing DFS analysis to assess frequency content,
   - Decimating the FEA data to a sampling rate that is similar to the
     test data,
   - Defining and assessing the FRF of a lowpass filter,
   - Applying the lowpass filter to both the test and FEA data.
   - Utilizing SRS analysis.
 o Utilizing Units with all the calculations, including using  
   k_unitsVariables to define some handy "units variables".
 o Using a tabbed figure to organize all the plots. 

Note: This example processes three data files. It will be most 
efficient to store these three data sets in a cell array and then loop
through the cell for each type of analysis process. We could also use
a structure for this, but this complicates the looping slightly and
is thus not used here. 

Tips:
 o Consider running the example "By Section". It will help you walk 
 through the example and its results.
%}

%% 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 inherited during
   certain math operations.
%}

k_cleanup();
adv = k_adv();
k_unitsPreferenceActivate('ShockVib_mm_ms_kG'); 
k_unitsVariables('mm, ms, N, G, kG, kHz');
k_set('inheritPropsAutoWarning', 'off');

% Define a tabbed figure to hold all the results
figH = figure('position', k_figZoom(2,2), ...
    'Name', 'MoonPenetrator');

titles = {
    'Raw Data'
    'Trimmed, Extended and Regularized'
    'DFS Assessment & Decimation'
    'Upsampling to Improve "Look of Data"'
    'Additional Lowpass Filtering'
    'Integrating for Velocity'
    'Shock Response Analysis via PVSS'
    'Final Summary of Results'    
    };

tabStuff = k_figTabsCreate(titles, figH);

%% Reading in the text files and displaying raw data.
%{
There are three text files for this example, each one is a simple 
file with a single column header row that also includes units. 

Additional Details: 
The commands below do the following:
 o Specify the data folder which contains the example's data and the
   filenames to use.
 o Reads the text data files in via k_readText.   
%}

dataDir = '..\data';

% Three files
fileNames = {
     'moon_testDataAll.csv'
     'moon_fea1.csv'
     'moon_fea2.csv'
     };
nsets = numel(fileNames);         

% Read in the data files
% Note: all three files are similar enough that we can use the same
% settings for the k_readText function on all the files 
delim = ',';
startStr = [];
rs1 = 0;
colHeadStr = 'time';
rs2 = 0;
endStr = [];
rs3 = [];

raw = cell(nsets, 1);

for z = 1 : nsets
    fName = fullfile(dataDir, fileNames{z});  % Path + file name
    raw{z} = k_readText(fName,delim,startStr,rs1,colHeadStr,rs2,...
        endStr,rs3);
end

%{
To get a quick view of the data, issue commands like this in the command
window:
   raw
   raw{1}
   k_varViewer(raw{1})
%}

% Display the results in Raw Data tab
currentTabH = tabStuff.TabsH(1);

% Show the data on the left side of the tab and plots on the right side.

dispBrief_orig = k_get('dispBrief');
k_set('dispBrief', 4);

ADVs = {
    'parent', currentTabH, ...
    'box', 'on'
    };

k_displayOnFigure(0.25, 5/6, raw{1}, ADVs);
k_displayOnFigure(0.25, 0.5, raw{2}, ADVs);
k_displayOnFigure(0.25, 1/6, raw{3}, ADVs);

k_set('dispBrief', dispBrief_orig);

for z = 1 : nsets
    subplot(3,2,2*z, 'parent', currentTabH)
    k_plot(raw{z}, [], ...
        'title', 'description') 
end

k_figTabsDisplay(currentTabH);

msgTxt = {
    'Notice the following:'
    ' o Test data has much longer time duration than FEA.'
    ' o FEA is much noisier and thus more freq content.'
    ' o FEA accel is in mm/s^2, whereas test is in G'
    };
uiwait(msgbox(msgTxt, 'Things to notice'));      

%% Trimming all data to common time duration and convert units
%{
This section trims all the data to a common time duration based on the
shortest duration among the three data sets. We also regularize each
data set to ensure that each has a constant sampling rate so that DSP
functions can be applied properly. Lastly, because the data starts 
abruptly, especially the experimental data, we pre-pad the beginning
of the data with zero amplitude.

Additional Details: 
 o The first thing done is to initialize the variable "endTimes" as a 
   k_units variable that will hold the max times of the three
   data sets. The "zeros(nsets,1)" produces a 3x1 array of zeros and
   then this is multiplied by the first value of time from raw{1},
   which is a simple "trick" to get a 3x1 array of 0 with the correct
   time units.
   - A similar trick is done for "fs" which will hold the sampling
     frequency of each of the data sets.
%}

% Initialize a couple of k_units variables 
endTimes = zeros(nsets,1)*raw{1}(1,'time');

fs = zeros(nsets,1)/raw{1}(1,'time');
fs = fs.convert('kHz');
fs.Props.Description = '';
fs.Props.ColNames = 'fs';
fs.Props.RowNames = fileNames;

for z = 1 : nsets
    endTimes(z,1) = max(raw{z}('time'));
end

% The maxCommonTime is the minimum of the endTimes
maxCommonTime = min(endTimes);

% Trim all data to this value of time. 
trimmed = cell(nsets,1);

% This prepad factor is set to help minimize filter start-up issues at
% at the beginning of the data due to its abrupt start. 
prePadFactor = 0.2; 

for z = 1 : nsets
    tmp = raw{z};
    idx = k_trimIndices(tmp('time'), [], maxCommonTime);
    tmp = tmp(idx,:);
    tmp = tmp.convert(); % Convert to preferred units
    tmp = k_regularize(tmp);
    tmp = k_extend(tmp, [], prePadFactor, 'zero', [], []);  
    fs(z,1) = k_samplingFreq(tmp('time'));
    trimmed{z} = tmp;
end

% Plot all of the data sets
currentTabH = tabStuff.TabsH(2);

pltTmp = cell(nsets,1);
for z = 1 : nsets
    subplot(3,1,z, 'parent', currentTabH)
    pltTmp{z} = k_plot(trimmed{z}, [], ...
        'title', 'description', ...
        'grid', 'both');
    
    txt = ['fs = ', num2str(fs(z,1))];
    k_displayOnFigure(0.8, 0.15, txt, ...
        'box', 'on');
end

k_plotMatchLimits(pltTmp, 'x');

k_figTabsDisplay(currentTabH);

msgTxt = {
    'Note the following:'
    ' o Magnitude of accel for FEA much larger than test '
    ' o FEA has much greater noise compared to test.'
    };
uiwait(msgbox(msgTxt, 'Things to notice'));  


%% Fourier Analysis and Decimation
%{
This section uses DFS (Discrete Fourier Spectrum) analysis and
decimation to get all the data sets to a bandwidth where a fair 
comparison can occur. 

Additional Details: 
o From the previous plots we see the max sampling frequency of the
  physical testing is 10*kHz. This means its maximum possible 
  frequency content is 5*kHz (its Nyquist frequency). This sampling
  is rather coarse for such a transient event, but we have been 
  assured that the physical test was recorded with a DAQ that has
  an antialiasing filter and that the data was further decimated
  using an antialiasing filter. 
  The FEA models have roughly 1000 times that frequency potential 
  due to their high sampling freq. 

o Below we compute a DFS of the data as-is and after it has been 
  decimated to a common sampling rate for all 3 datasets.  
%}

% Decimate FEA (the test is not being further decimated)
dec = cell(nsets, 1);
dec(1) = trimmed(1); % The test data is not additionally decimated
fsDec = fs; % Initialize fsDec

for z = 2 : nsets
    tmp = trimmed{z};
    decFactor = fs(z,1)/fs(1,1); % note: fs(1,1) = 10*kHz
    
    dec{z} = k_decimate(tmp, [], decFactor, ...
        adv.k_decimate.assumeStart.zero, ...
        adv.k_decimate.decimateFactorMax, 1200); % Extend the decimateFactorMax slightly
    dec{z}.Props.Description = tmp.Props.Description;
    
    fsDec(z,1) = k_samplingFreq(dec{z}('time'));
end


% Compute DFS BEFORE and after decimating FEA
dfsBefore = cell(nsets,1);
dfsAfter = cell(nsets,1);

for z = 1 : nsets    
    % Before decimation
    tmp = trimmed{z};
    dfsBefore{z} = k_DFS(tmp, [], 'hann', 'Mag');
    dfsBefore{z}.Props.Description = tmp.Props.Description;
    
    % After decimation
    tmp = dec{z};
    dfsAfter{z} = k_DFS(tmp, [], 'hann', 'Mag');
    dfsAfter{z}.Props.Description = tmp.Props.Description; 
end

% Plot the DFS results
currentTabH = tabStuff.TabsH(3);

subplot(2,2,1, 'parent', currentTabH);
k_plot(dec, [], ...
    'legendLocation', 'SouthEast', ...
    'title', {'user', 'Time-Domain of All 3 with Common fs'}) 


plotADVS = {'lineWidth', [2,0.5, 0.5], ...
            'legendLocation', 'NorthEast'};

subplot(2,2,2)
k_plot(dfsBefore, [], ...
    plotADVS, ...
    'scale', 'logLog', ...
    'xLim', [0.1,100]*kHz, ...
    'legendLoc', 'SouthEast', ...
    'title', {'user', 'DFS of test and FEA BEFORE Decimating FEA'}) 

subplot(2,2,4)
k_plot(dfsAfter, [], ...
    plotADVS, ...
    'title', {'user', 'DFS of test and FEA AFTER Decimating FEA'}) 

msgTxt = {
    'Notes: '
    ' '
    '1) Both the time-domain and DFS plots show that the fea1 model is ' 
    '   the better match to testData.'
    ' '
    '2) The decimated FEA and the testData look a bit coarsely sampled '
    '   at 10*kHz sampling frequency. '
    };

k_displayOnFigure(0.1, 0.4, msgTxt, ...
    'parent', currentTabH, ...
    'box', 'on', ...
    'alignHor', 'left', ...
    'alignVert', 'top')

k_figTabsDisplay(currentTabH);

%% Upsample the results to improve the "look" of the data
%{
Due to the frequency content having noticeable magnitude modestly near
the Nyquist frequency, the time-domain signals look a bit "discrete".
To address this, we upsample the signals. This does not add frequency
content, it just represents the data with more points using a 
DSP-based upsampling filter. We also need to upsample for the 
SRS analysis that will be performed later. 
%}

Q = 50; % upsampling factor.

up = cell(nsets,1);
fsUp = fsDec;
for z = 1 : nsets
    up{z} = k_upsample(dec{z}, [], Q);
    fsUp(z,1) = k_samplingFreq(up{z}('time'));    
end

% Plot the time-domain results due to upsampling
currentTabH = tabStuff.TabsH(4);

k_plot(up, [], ...
    'parent', currentTabH, ...
    'legendLocation', 'SouthEast', ...
    'grid', 'both', ...
    'title', {'user', 'Time-domain of all 3 curves after upsampling'})

txt = 'Model fea1 is clearly the best match to the physical test';
k_displayOnFigure(0.5, 0.2, txt, ...
    'box', 'on', ...
    'fontSize', 14, ...
    'fontWeight', 'bold');

k_figTabsDisplay(currentTabH);

%% Apply additional LP filtering to further compare
%{
Below we apply an additional lowpass filter to compare the FEA1 
model and the physical test. It is noted that the calculation of the
FRF is optional, but it aids in understanding what the filter is
doing. 

Note, since all the upsampled curves have the same sampling frequency, we
can compute a single set of filter coefficients using fsUp(1,1).
%}

% Lowpass filter cutoff frequency
fc = 1.5*kHz;

% Calculate filter coefficients 
c = k_cIIR(fsUp(1,1), fc, 'lp', 'but');
nPasses = 2; % Bi-directional filter

% FRF calcs
fVec = k_vectorByPoints(0, 1, 100)*5*kHz; 
frf = k_FRF(fsUp(1,1), fVec, c, nPasses, 'Mag');   
frf.Props.Description = ['FRF for fc = ', num2str(fc)];

nResults = 2;
results = cell(nResults,1);
for i = 1 : nResults        
    results{i} = k_filter(up{i}, [], c, nPasses, ...
        'assumeStart', 'zero'); 
end


% Plot FRF and Filtered results
currentTabH = tabStuff.TabsH(5);

subplot(2,2,1, 'parent', currentTabH)
k_plot(dfsAfter(1:2), [], ...
    plotADVS, ...
    'xGuideVal', fc, ...
    'title', {'user', 'DFS of test and FEA1 AFTER Decimating FEA'}) 

subplot(2,2,2)
k_plot(up(1:2), [], ...
    'markerSymbols', 'list', ...
    'markerOccurrence', {'total', 7}, ...
    'legendLoc', 'SouthEast', ...
    'title', {'user', 'Upsampled Results before LP filtering'}) 

subplot(2,2,3)
k_plot(fVec, frf, ...
    'xLim', [0,1]*5*kHz, ...
    'xGuideVal', fc, ...
    'legendLocation', 'NorthEast', ...
    'title', 'yDescription') 

subplot(2,2,4)
titleTxt = k_stringCreate({'Filtered Result for fc = ', fc});
k_plot(results, [], ...
    'markerSymbols', 'list', ...
    'markerOccurrence', {'total', 7}, ...
    'legendLoc', 'SouthEast', ...
    'title', {'user', titleTxt}) 

k_figTabsDisplay(currentTabH);

%% Integrating Accel to assess velocity
%{
It is generally good practice to integrate shock acceleration data
to make sure the velocity data looks plausible.

In this example of a penetrator, the initial velocity was not supplied
so the calculation of velocity will be off by a constant value, that of
the initial velocity.
%}

vel_up = cell(nResults,1);
vel_results = cell(nResults,1);

for z = 1 : nResults
    vel_up{z} = k_integrate(up{z});
    vel_up{z}.Props.ColNames{2} = 'velocity';
    
    vel_results{z} = k_integrate(results{z});    
    vel_results{z}.Props.ColNames{2} = 'velocity';    
end

% Plot results
currentTabH = tabStuff.TabsH(6);

subplot(2,1,1, 'parent', currentTabH)
k_plot(vel_up, [], ...
    'legendLocation', 'NorthEast', ...
    'grid', 'both', ...
    'title', {'user', 'Velocities of Upsampled Accelerations'})

subplot(2,1,2, 'parent', currentTabH)
titleText = ['Velocities of LP Filtered Accelerations w/ fc = ', num2str(fc)];
k_plot(vel_results, [], ...
    'legendLocation', 'NorthEast', ...
    'grid', 'both', ...
    'title', {'user', titleText})

k_figTabsDisplay(currentTabH);

%% Shock Response Analysis
%{
A good way to compare transient shock data is via SRS analysis. 
Below we compute a PVSS (Pseudo Velocity Shock Spectrum) to compare
the Physical Test to the FEA.
%}

damp = 0.05; % Assumped damping of the payload. 
pvss = cell(nResults,1);

for z = 1 : nResults
    pvss{z} = k_SRS(up{z}, [], '1*Hz', 'fnMax', damp, 'pvss', ...
        'modifyColNames', 'replace'); 
end

% Plot results
currentTabH = tabStuff.TabsH(7);
k_plot(pvss, [], ...
    'parent', currentTabH, ...
    'lineWidth', 2, ...
    'scale', 'fourCP', ...
    'LegendLoc', 'NorthEast', ...
    'title', 'yDescr', ...
    adv.k_plot.descripRowForTitle.last, ...
    'yLim', '[-inf, 1e3]*mm/ms')

k_figTabsDisplay(currentTabH);

%% Making a final summary figure
%{
In this section we make a single final summary figure that nicely
summarizes the results we achieved.
%}

currentTabH = tabStuff.TabsH(end);

subplot(3,2,1, 'parent', currentTabH)
k_plot(raw{1}, [], ...
    'grid', 'both', ...
    'units', 'ms, kG', ...
    'title', {'user', 'Test - Raw Data'})

subplot(3,2,2)
k_plot(raw{2}, [], ...
    'grid', 'both', ...    
    'yLim', [-50, 50]*kG, ...
    'units', 'ms, kG', ...    
    'title', {'user', 'FEA1 - Raw Data'})

subplot(3,2,3)
k_plot(dfsBefore(1:2), [], ...
    'legendText', {'user', {'Test', 'FEA1'}}, ...
    'grid', 'both', ...    
    'xLim', [0, 6]*kHz, ...    
    'legendLocation', 'NorthEast', ...
    'title', {'user', 'Fourier Spectrum'})    

subplot(3,2,4)
k_plot(fVec, frf, ...
    'yLim', [0,1.1], ...
    'xGuideVal', fc, ...
    'xGuideWidth', 2, ...
    'title', {'user', 'Filter''s FRF'})   

k_displayOnFigure(fc*1.1, 0.8, fc, ...
    adv.k_displayOnFigure.locationUnits.mixedDataAndNorm, ...
    'removeVarName', 'off', ...
    'alignHor', 'left')

subplot(3,2,5)
k_plot(results{1}, [], ...
    'grid', 'both', ...    
    'yLim', [-6,2]*kG, ...
    'title', {'user', 'Test - After Filtering'})

subplot(3,2,6)
k_plot(results{2}, [], ...
    'grid', 'both', ...    
    'yLim', [-6,2]*kG, ...
    'title', {'user', 'FEA1 - After Filtering'})

k_figTabsDisplay(currentTabH);

%% Options to save figure windows

k_exampleSetup('figSavingOptions');
