function [newburst,qc]=CookRawADV(burst,nsd,ndt,min_amp,min_corr,b_title,wantplot,CHOICE); 
% Identify and remove suspicious points from ADV data
% [newburst,qc]=CookRawADV(burst,nsd,ndt,min_amp,min_corr,wantplot,CHOICE);
% 
% This code enables the user to remove spurious data points from ADV
% time series.  5 criteia are employed to flag bad data and the user
% will be notified of the bad points via a colour-coded plot.
%
% The input arguments are as follows with some recommended initital values:
% nsd=2.8;        % num std. deviations that defines outliers
% ndt=6;          % cutoff frequency
% min_amp=120;    % Amplitude Threshold (Worst average is ~127)
% min_corr=60;    % Correlation Threshold
% b_title='string'; % The Burst Number ID comes from the Header in ADV_Bbq.m
% wantplot=1;     % If you wish to view diagnostic plots, wantplot=1
% CHOICE=[1 1 1 1 1]; % For each test 1 through 5, 1=yes >>DEFAULT IS ALL YES<<
% 
% CHOICE LEGEND >>>
% #1 = signal strength test => green dot
% #2 = correlation test => magenta dot
% #3 = rapid shift in slope test => cyan circle
% #4 = slope deviates strongly from local slope => green cross
% #5 = fails low pass filter => red x
%
% Originally Written by Chris Sherwood, USGS
%
% Restructured and revised 9/15/03 by Jason M. Fox for Paul Hill

% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% Establish A Stipulation whereby the CHOICE options can be modified
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
burst_ok=[0];

while burst_ok==0
    
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% Stipulation is at WORK!
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    % Before Cooking the Raw ADV, Initialization must take place
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    
    % Create the Output Data Structure
    newburst=burst;
    
    % number of points
    n = length(burst.Vx{1});
    fs = burst.samprate;
    delta_t = 1/fs;
    cutoff_f = 0.1*fs;
    
    % Mirror Original Velocity data
    ur = burst.Vx{1};
    vr = burst.Vy{1};
    wr = burst.Vz{1};
    
    % Build Arrays for Replacement values (low-pass filtered somthin')
    urr = burst.Vx{1};
    vrr = burst.Vy{1};
    wrr = burst.Vz{1};
    
    % Build an array for storing bad data flags
    b = zeros(n,6);
    
    % Has the user input selected CHOICE?  If no, set default CHOICE
    if exist('CHOICE')==0
        CHOICE=[1 1 1 1 1];
    end
    
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    % Initialization Complete
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    % Compile Quality Flags
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    
    %% bad test 1 - flag poor signal strength ratio (threshold test)
    if CHOICE(1)==1
        b(:,1) = (burst.amp{1}(:,1)<min_amp)+(burst.amp{1}(:,2)<min_amp)+(burst.amp{1}(:,3)<min_amp);
    end
    
    %% bad test 2 - flag poor correlation (threshold test)
    if CHOICE(2)==1
        b(:,2) = (burst.corr{1}(:,1)<min_corr)+(burst.corr{1}(:,2)<min_corr)+(burst.corr{1}(:,3)<min_corr);
    end
    
    %% bad test 3 - flag big changes in u or v
    % base slope test with combined absolute diff of detrended u and v time series
    if CHOICE(3)==1
        DU = abs(diff(detrend(ur)));
        DV = abs(diff(detrend(vr)));
        biguslp = DU>4*std(DU);
        bigvslp = DV>4*std(DV);
        
        % pairs of big DX indicate an outlier.
        % Flag 'em!
        b(:,3) = ( ([0; biguslp] == [biguslp; 0]) & [biguslp;0]~=0 ) + ...
            ( ([0; bigvslp] == [bigvslp; 0]) & [bigvslp;0]~=0 );
        
        %% bad test 4 - flag when slope differs from local slope
        DUf = lpfilt(DU,delta_t,cutoff_f);
        DVf = lpfilt(DV,delta_t,cutoff_f);
        DUr = DU-DUf;
        DVr = DV-DVf;
        bigdus = abs(DUr)>3*std(DUr);
        bigdvs = abs(DVr)>3*std(DVr);
    end
    
    % pairs of outlier DX indicate an outlier point
    %Flag 'em!
    if CHOICE(4)==1
        b(:,4) = ( ([0; bigdus] == [bigdus; 0]) & [bigdus;0]~=0 ) + ...
            ( ([0; bigdvs] == [bigdvs; 0]) & [bigdvs;0]~=0 );
    end
    
    %% bad test 5 - sticks out from smoothed time series
    % low-pass filter
    if CHOICE(5)==1
        urf=lpfilt(ur,delta_t,cutoff_f);
        vrf=lpfilt(vr,delta_t,cutoff_f);
        wrf=lpfilt(wr,delta_t,cutoff_f);
        uresid = ur-urf;
        vresid = vr-vrf;
        
        % Keep in mind other flags could be used to find replacement points.
        % The low-pass filter works well for the time being.
        bad = ( (abs(uresid) > nsd*std(uresid)) | ( abs(vresid) > nsd*std(vresid) ));
        b(:,5)=bad;
        rep = find(bad);
        nrep = find(~bad);
        
        % generate a revised velocity profile with replacement from the lpfilt profile
        if(any(rep)),
            urr(rep)=urf(rep);
            vrr(rep)=vrf(rep);
            wrr(rep)=wrf(rep);
        end
    end
    %% bad test summation 6 - note flagged points for tests 1 through 5 in the end column
    b(:,6) = sum( b(:,1:5)')';
    
    % Now that we know all the good guys from the first 4 test, let's establish
    % a profile for a trace sans flagged dudes
    bee_in_bonnet=find(b(:,6)==0);
    
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    % Quality Flags Complete
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    % Create a revised Series through Linear Interpolation
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
            
    % Define the Time Frame for the series    
    bt = (0:delta_t:n*delta_t-delta_t)';
    
    new_u=interp1(bt(bee_in_bonnet),ur(bee_in_bonnet),bt);
    new_v=interp1(bt(bee_in_bonnet),vr(bee_in_bonnet),bt);
    new_w=interp1(bt(bee_in_bonnet),wr(bee_in_bonnet),bt);
    
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    % Build a Quality Report
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    
    % recalculate mean and variance
    qc.badpts = sum(b);
    if exist('rep')==0
        qc.nfixed = 0;
    else
        qc.nfixed = length(rep);
    end
    qc.old_mean  = mean([ ur, vr, wr]);
    qc.new_mean_linear  = mean([ new_u, new_v, new_w ]);
    qc.new_mean_lowpass  = mean([ urr, vrr, wrr ]);
    qc.old_var   = std([ ur, vr, wr ]).^2;
    qc.new_var_linear   = std([ new_u, new_v, new_w ]).^2;
    qc.new_var_lowpass   = std([ urr, vrr, wrr ]).^2;
    qc.old_max = max([ ur, vr, wr ]);
    qc.new_max_linear = max([ new_u, new_v, new_w ]);
    qc.new_max_lowpass = max([ urr, vrr, wrr ]);
    qc.old_min = min([ ur, vr, wr ]);
    qc.new_min_linear = min([ new_u, new_v, new_w ]);
    qc.new_min_lowpass = min([ urr, vrr, wrr ]);
    
    % Post the Report in the Command Window
    qc
    
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    % Quality Control Report Complete
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    % Prepare and Post Plots for Review
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    
    if wantplot==1 
        
        % Prepare Detrended data for the plots
        nd(:,1) = detrend(ur)./std(detrend(ur));
        nd(:,2) = detrend(vr)./std(detrend(vr));
        ndr(:,1) = detrend(urr)./std(detrend(ur));
        ndr(:,2) = detrend(vrr)./std(detrend(vr));
        ndf(:,1) = lpfilt(nd(:,1),delta_t,cutoff_f);
        ndf(:,2) = lpfilt(nd(:,2),delta_t,cutoff_f);
        
        % Set the appearance of the symbols for each quality test
        sym = ['.g';'.m';'oc';'+g';'xr'];
        
        % Set the number of time steps per plot
        npts=200; % at 10 Hz, this gives 20 s.
        
        
        for ich=1:n/npts;
            % Plot U V (Green Solid LINE), smoothed detrended U V (Blue Dashed LINE),
            % detrended UV points (black POINTS)
            % and appropriately flagged points (ASSORTED SYMBOLS)
            clf
            ie = npts*ich;
            is = ie-npts+1;
            
            % Make suntin' outta nuttin'
            nuttin=[NaN];
            
            % Plot the Full burst as the top slide            
            subplot(311)
            [AX,H1,H2]=plotyy(bt,ur(:,1),nuttin,nuttin);
            hold on;
            % Title
            title(b_title,'Fontsize',12);
			% Set Axes Ranges
			set(AX(1),'Xlim',[bt(1) bt(length(bt))]);
			set(AX(1),'Ylim',[-40 40]);
			set(AX(2),'Xlim',[bt(1) bt(length(bt))]);
			set(AX(2),'Ylim',[-40 40]);
			% Place the proper tick marks on Plot
			set(AX(1),'Ytick',[-40 -30 -20 -10 0 10 20 30 40],'YtickLabel','-40| |-20| |0| |20| |40');
			set(AX(2),'Ytick',[-40 -30 -20 -10 0 10 20 30 40],'YtickLabel','-40| |-20| |0| |20| |40');
			% Y-axes labels
			set(get(AX(1),'Ylabel'),'String','U Velocity (cm/s)','fontsize',9);
			set(get(AX(2),'Ylabel'),'String','V Velocity (cm/s)','fontsize',9);
            % Now mark the section of the burst in view in the panels below 
            plot(bt,vr(:,1),'-g');
            plot(bt(is:ie),vr(is:ie,1),'-k',bt(is:ie),ur(is:ie,1),'-r');
            
            % U velocity
            subplot(312)
            % Start with filtered detrended U on the left
            % and U velocity on the right
            [AX,H1,H2]=plotyy(bt(is:ie),ur(is:ie,1),bt(is:ie),ndf(is:ie,1));
            hold on;
			% Set Axes Ranges
			set(AX(2),'Xlim',[bt(is) bt(ie)]);
			set(AX(2),'Ylim',[-4 4]);
			set(AX(1),'Xlim',[bt(is) bt(ie)]);
			set(AX(1),'Ylim',[-40 40]);
			% Place the proper tick marks on Plot
			set(AX(2),'Ytick',[-4 -3 -2 -1 0 1 2 3 4],'YtickLabel','-4| |-2| |0| |2| |4');
			set(AX(1),'Ytick',[-40 -30 -20 -10 0 10 20 30 40],'YtickLabel','-40| |-20| |0| |20| |40');
			% Y-axes labels
			set(get(AX(2),'Ylabel'),'String','Std. Dev. about Filtered Detrended Ur','fontsize',9);
			set(get(AX(1),'Ylabel'),'String','U Velocity (cm/s)','fontsize',9);
            % Now add the series points and new interpolated U
            plot(bt(is:ie),ur(is:ie,1),'.k')
            plot(bt,new_u,'-r')
            % Finally, add the flags
            for ibad=1:5,
                blist = find( b(:,ibad) );
                if(any(blist))
                    plot(bt(blist),ur(blist,1),sym(ibad,:));
                end
            end
            
            % V velocity
            subplot(313)
            [AX,H1,H2]=plotyy(bt(is:ie),vr(is:ie,1),bt(is:ie),ndf(is:ie,2));
            hold on;
            % Set Axes Ranges
			set(AX(2),'Xlim',[bt(is) bt(ie)]);
			set(AX(2),'Ylim',[-4 4]);
			set(AX(1),'Xlim',[bt(is) bt(ie)]);
			set(AX(1),'Ylim',[-40 40]);
			% Place the proper tick marks on Plot
			set(AX(2),'Ytick',[-4 -3 -2 -1 0 1 2 3 4],'YtickLabel','-4| |-2| |0| |2| |4');
            set(AX(1),'Ytick',[-40 -30 -20 -10 0 10 20 30 40],'YtickLabel','-40| |-20| |0| |20| |40');
			% Y-axes labels
			set(get(AX(2),'Ylabel'),'String','Std. Dev. about Detrended Vr','fontsize',9);
			set(get(AX(1),'Ylabel'),'String','V Velocity (cm/s)','fontsize',9);
            % X-Axis Label
            set(get(AX(1),'Xlabel'),'String','Time (s)','fontsize',9);
            % Now add the series points and new interpolated V
            plot(bt(is:ie),vr(is:ie,1),'.k')
            plot(bt,new_v,'-r')
            % Finally, add the flags
            for ibad=1:5,
                blist = find( b(:,ibad) );
                if(any(blist))
                    plot(bt(blist),vr(blist,1),sym(ibad,:));
                end
            end
            
            
            drawnow
            shg
            
            % Allow the user to take a pause and look closely
            pause
        end
    end
    
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    % Plots Complete
    % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    
    % Ask the user if he/she accepts the revised set or wishes to revise the
    % CHOICES options
    
    %% >>>>>>>>>>>>>BEGIN QUESTION
    ButtonName=questdlg('Is the revised data set Acceptable or do you wish to revise your test options?  NOTE: If you choose to revise your options, the command prompt will await your input.', ...
        'Is the ADV data cooked to your liking?', ...
        'Accept','Revise CHOICE options','Accept');
    
    switch ButtonName,
        case 'Accept', 
            burst_ok=1;
        case 'Revise CHOICE options',
            burst_ok=0;
    end
    %% >>>>>>>>>>>>>>END QUESTION
    
    if burst_ok==1
        break;
    end
    
    %% Now Re-define the Choices Options for the upcoming recalculation
    CHOICE = input('Enter the revised CHOICE options set (include sq. brackets [x x x x x]) : ');
end

% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% Prepare the Revised Data Set for Export
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

newburst.Vx{:}=new_u;
newburst.Vy{:}=new_v;
newburst.Vz{:}=new_w;

% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% Done
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>