vectrinoDataScrub.m

function [U,badInds]=vectrinoDataScrub(Data,Config,opts)
% Post processing code for Vectrino ADV velocity profiles
% Input, Data and Config are the structs generated by conversion of .ntk
% vectrino files to MATLAB format using the Vectrino software.  Velocity
% can be measured in either beam (preffered) or xyz coordinates
% NOTE: Currently this is not setup to proccess Vectrino data in profiling
% mode (more than 1 interrogation cell). This functionality may be added in
% the future

% opts =
%   0: (default) replaces bad points with NaN
%   1: Removes bad points
%   2: Replaces bad points with polynominal interpolation via neighbors

% This script removes points with "bad" correlation, either every point
% with correlation below 70% or the 10% of points with the worst
% correlation, whichever results in the removal of the fewest points

% In addition, despiking is performed on the beam velocities.  However,
% a high pass filter applied prior to despiking, reducing variance due to
% slow velocity oscillations and incresing the effectiveness of the
% despiking algorithm.

if ~exist('opts','var') || isempty(opts)
    opts=0;
end
%% If needed, convert back to beam coordinates

% Get the tranform matrix
T= reshape(Config.ProbeCalibration_calibrationMatrix,4,4)';

% Get the velocity, transform if needed
if isfield(Data,'Profiles_VelX')
    Uin=[Data.Profiles_VelX';Data.Profiles_VelY';...
        Data.Profiles_VelZ1';Data.Profiles_VelZ2'];
    Ubeam=(T\Uin)';
else
    Ubeam=[Data.Profiles_VelBeam1,Data.Profiles_VelBeam2,...
        Data.Profiles_VelBeam3,Data.Profiles_VelBeam4];
end

%% Find bad points based on correlation
% NOTE: If a point is bad in one beam, all beams are considered bad.
cor=[Data.Profiles_CorBeam1, Data.Profiles_CorBeam2, Data.Profiles_CorBeam3, Data.Profiles_CorBeam4];
corLogic=ones(size(cor));
corLogic(cor<70)=0;
corLogic2=prod(corLogic,2);
badCor=zeros(length(corLogic2),1);
if sum(corLogic2==0)<0.1*length(corLogic2)
    % Then chuck the ones less than 70
    badCor(corLogic2==0,:)=1;
else
    % Chuck the worst 10%
    [~,inds]=sort(cor);
    inds=reshape(inds',numel(inds),1); 
    uniqueInds=unique(inds,'stable');
    chuckInds=uniqueInds(1:round(length(corLogic2)*0.1));
    badCor(chuckInds)=1;
end

%% Find bad points based on despike
sf=Config.sampleRate;
cutoffFreq=3; % High pass filter cuttoff in Hz. May need to adjust
badSpike=zeros(length(corLogic2),1);
dif = designfilt('highpassfir', 'FilterOrder', 20, 'CutoffFrequency', cutoffFreq, 'SampleRate', sf);
for i=1:4
    
    f=filtfilt(dif,detrend(double(Ubeam(:,i))));
    
%     [~,spikeInds]=func_despike_phasespace3d(f);
%     [~,spikeInds]=func_despike_phasespace3d(f, 9, 0);
    [~,spikeInds]=func_despike_phasespace3d(f, 9, 2);
      
    badSpike(spikeInds)=1;
end
    
%% Merge correlation and despike
badInds=zeros(size(badSpike));
badInds(badSpike | badCor)=1;
badInds=logical(badInds);
%% If needed, interpolate/remove points. Otherwise NaN
x=find(~badInds);
xq=find(badInds);
if opts==1
    Ubeam(badInds,:)=[];
elseif opts==0
    Ubeam(badInds,:)=NaN;
elseif opts==2
    for i=1:4
        Ubeam(xq,i)=interp1(x,Ubeam(x,i),xq,'PCHIP');
    end
end

%% Transform back to XYZ choords
U=(T*Ubeam')';
end