function [  ] = SOME_pitch_testing(wav,fs,ind,context)
%SOME_PITCH_TESTING Tests parameters for fxrapt, outlier removal &
% smoothing. Parameters for testing should be arrays (of doubles); the
% first value should be the default value (e.g., as given in the fxrapt
% documentation), for plotting purposes (or modify the plot title).
%%%
%   The backbone is Sam Tilsen's contour_processing.m, modified to take 
%   one wav file at a time. The processing functions are from
%   ST_pitch_process.m, in Bruce's pitch visualization GUI.
%%%
%   Input arguments: 
%       wav = the audio for pitch contour processing
%       fs = sampling rate of audio file
%       ind = the EZRA index of the wav file, for plotting
%       context = the word following SOME, for plotting

%% parameters
d = fdesign.highpass('N,Fp,Ast,Ap',3,70,3,0.8,fs);
Hd = design(d,'ellip'); % Starts FVTool to display the filter.
reset_voicebox_rapt;

%%% Semi-fixed fxrapt params (I don't expect them to make a big difference)
% voicebox('rapt_tframe', 0.0075);  % frame size for fxrapt (from ST); default 0.01
% voicebox('rapt_nfullag', 15);   % default is 7
% voicebox('rapt_tcorw',0.02);  % doesn't do anything? correlation window
% size; default 0.0075
% voicebox('rapt_vtrac',.9); % default 0.5 (delta amplitude modulated transition cost)
% voicebox('rapt_absnoise', 0.05);    % default 0
% voicebox('rapt_relnoise', 2);
% voicebox('rapt_signoise',.0005); % default .001 (ratio of peak signal rms to noise floor)

%%% Params for post-processing
par.max_F0_jump = 15;            % outlier removal param
par.sd_tol = 2.4;                       % outlier removal param
par.max_frames_interp = 4;           % interpolation param;; GUI fct doesn't use it
par.sm_bw = 5;        % parameter for smoothing
par.sm_params = [0.95 .995];  % smoothing params
%par.sm_params = [0.95 .9995];

%%% fxrapt params for testing; NOTE: 1st value is the default
par.doublec = [0.35 0.7 1 1.5];           % added by AC for halving/doubling, but it doesn't seem to be working
% par.doublec = [0.35 1.5 3 5];
par.vtranc=[0.005 0.01 0.2]; % fixed voice-state transition cost (default .005)
par.freqwt=[0.02 0.04 0.08 0.1]; % cost factor for F0 change (default 0.02; gets rid of spikes)
par.nfullag=[7 9 15]; % ST has it at 15.
par.tcorw=[.0075 .02]; % correlation window size. ST has it at .02
par.vtrac=[.5 .7 .9 1.3]; % delta amplitude modulated transition. ST has it at .9
par.absnoise=[0 .03 .05 .09]; % ST has it at .05
%par.relnoise=
par.signoise=[.001 .0005 .005 .01]; % ratio of peak signal rms to noise floor. ST has it at .0005

%%% remove_F0_outliers params for testing
% par.sd_tol = [2.0 2.2 2.4 2.6];
% par.max_F0_jump = [8 10 15 17];

%% atm I do not have gender encoded, so I'm hard-coding this
par.gend = 'M';
switch(par.gend)
    case 'F'
        voicebox('rapt_f0min',125);
        voicebox('rapt_f0max',400);
        
    case 'M'
        voicebox('rapt_f0min',75);
        voicebox('rapt_f0max',300); %initially 250, but in my data they seem to go up to 300
end

%% perform calculations

% I don't know what all of this wav pre-processing is for. Plotting the old
% wav (which came as input here) and the new wav (calculated in the next 4
% lines), I don't see any difference
wav = wav-mean(wav);
wav = filtfilt(Hd.sosMatrix,Hd.scaleValues,wav);
wav = wav-mean(wav);
wav = wav/max(abs(wav));
% this one though, it seems to be for plotting the waveform
twav = linspace(0,length(wav),length(wav))/fs;

%test_3_params('doublec', par.doublec, 'vtranc', par.vtranc, 'freqwt', par.freqwt);
test_3_params('vtrac', par.vtrac, 'absnoise', par.absnoise, 'signoise', par.signoise);
%test_3_params('doublec', par.doublec, 'nfullag', par.nfullag, 'tcorw', par.tcorw);
%test_remove_outliers(par.sd_tol, par.max_F0_jump);

%% FUNCTIONS START HERE
%% test 3 params
    function [y] = test_3_params(par1_name, par1, par2_name, par2, par3_name, par3)
 % tries all combinations of the 3 voicebox parameters. Displays and prints
 % a figure for each combination of par1 & par2, with multiple subplots for
 % par3. 
 % Input arguments:
 %      par1, par2, par3 = arrays of doubles. These are values you want to
 %      test for each parameter. The first value should be the default
 %      value for the plot titles to be accurate (or change plot title).
 %      par1_name, par2_name, par3_name = strings of the names of the
 %      voicebox parameters we're testing (without the "rapt_" part). They are used for plot titles,
 %      but also to update voicebox.
 
        rapt_par1 = ['rapt_',par1_name];
        rapt_par2 = ['rapt_',par2_name];
        rapt_par3 = ['rapt_',par3_name];
        
        % used to tell Matlab how many subplots there will be in 1 figure. 1 row
        % for the waveform; each contour will have it outlier-removed version
        % underneath.
        N_subplots = length(par3) + 1;
    
    
        for i=1:length(par1)
            voicebox(rapt_par1,par1(i));
            
            for j=1:length(par2)
                voicebox(rapt_par2, par2(j));
                h = figure; set(h, 'Units','inches','position', [5 11 17 15], 'PaperPositionMode','auto');
                
                % tight_subplot is not a Matlab function, but it allows me to
                % manipulate gaps and margins much easier. With Matlab's subplot
                % function, I couldn't make the contour subplots bigger without
                % having them overlap.
                ax = tight_subplot(N_subplots,2,[.05 .03],[.03 .01],[.03 .01]);
                
                % waveform. I'm making this small, since it's just for reference.
                axes(ax(1)); plot(twav,wav,'color',[.5 .5 .5]);
                p = get(ax(1),'pos'); p(4) = p(4) - 0.1; set(ax(1),'pos',p);
                axes(ax(2)); plot(twav,wav,'color',[.5 .5 .5]);
                p = get(ax(2),'pos'); p(4) = p(4) - 0.1; set(ax(2),'pos',p);
                
                for k=1:length(par3)
                    
                    voicebox(rapt_par3,par3(k));
                    
                    [ffr,t] = fxrapt(wav, fs, 'u');
                    ffrt = mean(t(:,1:2),2)/fs;
                    ffr_out = remove_F0_outliers(ffr,par.sd_tol,par.max_F0_jump);
                    
                    % figure out where to put the plots. Outlier plots underneath
                    % raw contours.
                    if rem(k,2)==1
                        pos_raw = 2*k+1;
                        pos_out = 2*k+3;
                    else
                        pos_raw = 2*k;
                        pos_out = 2*k+2;
                    end
                    
                    axes(ax(pos_raw));
                    plot(ffrt, ffr, 'o-'); hold on;
                    plot(ffrt(isnan(ffr)), min(ffr), 'ro');
                    title(['Raw, ', par3_name, ' = ', num2str(par3(k))],'Fontsize',12);
                    
                    axes(ax(pos_out));
                    plot(ffrt, ffr_out, 'go-'); hold on;
                    plot(ffrt(isnan(ffr_out)), min(ffr_out), 'ro');
                    title('Outliers removed, sd = 2.4, max-F0-jump = 15','Fontsize',12);
                end
                
                linkaxes(ax,'x');
                linkaxes(ax(3:length(ax)),'y');
                
                % Figure title, centered
                set(gcf,'NextPlot','add');
                axes;
                if i == 1
                    defi = 'DFLT';
                else
                    defi = '';
                end
                
                if j ==1
                    defj = 'DFLT';
                else
                    defj = '';
                end
                
                h2 = title(['SOME  ', upper(context), ' ', int2str(ind), ' raw (', par1_name, ' = ', defi, num2str(par1(i)), ...
                    ', ', par2_name, ' =  ', defj, num2str(par2(j)), ') and outlier removed (sd = 2.4, max-F0-jump = 15).'], 'Fontsize', 18);
                set(gca,'Visible','off');
                set(h2,'Visible','on');
                
                % sound(wav,fs);
                
                formatOut = 'mm.dd HHMM';
                timestmp = datestr(now,formatOut);
                filename = [timestmp, ' - ', int2str(ind), ' raw & outlier removed, ', par1_name, ' = ', num2str(par1(i)), ', ', par2_name, ' =  ', ...
                    num2str(par2(j)), ', ', par3_name, ' = ', regexprep(num2str(par3),' *',','),'.eps'];
                print(h, '-depsc', filename);
                
                % pause;
                close(h);
            end
        end
    end
%% test remove outliers
    function [y] = test_remove_outliers(sd_tol,max_F0_jump)
    
    N_subplots = floor(length(max_F0_jump)/2) + 2;
    
    % current best fxrapt parameter values
    % voicebox('rapt_tframe', 0.0075);  % frame size for fxrapt (from ST); default 0.01
    % voicebox('rapt_nfullag', 15);   % default is 7
    % voicebox('rapt_tcorw',0.02);  % doesn't do anything? correlation window
    % size; default 0.0075
    % voicebox('rapt_vtrac',.9); % default 0.5 (delta amplitude modulated transition cost)
    % voicebox('rapt_absnoise', 0.05);    % default 0
    % voicebox('rapt_relnoise', 2);
    % voicebox('rapt_signoise',.0005); % default .001 (ratio of peak signal rms to noise floor)
    voicebox('rapt_doublec',0.7);
    voicebox('rapt_freqwt',0.1);
    [ffr,t] = fxrapt(wav, fs, 'u');
    ffrt = mean(t(:,1:2),2)/fs;
    
    for k=1:length(sd_tol)
    
        h = figure; set(h, 'Units','inches','position', [5 11 17 15], 'PaperPositionMode','auto');
        
        % tight_subplot is not a Matlab function, but it allows me to
        % manipulate gaps and margins much easier. With Matlab's subplot
        % function, I couldn't make the contour subplots bigger without
        % having them overlap.
        ax = tight_subplot(N_subplots,2,[.05 .03],[.03 .01],[.03 .01]);
        
        % waveform. I'm making this small, since it's just for reference.
        axes(ax(1)); plot(twav,wav,'color',[.5 .5 .5]);
        p = get(ax(1),'pos'); p(4) = p(4) - 0.1; set(ax(1),'pos',p);
        axes(ax(2)); plot(twav,wav,'color',[.5 .5 .5]);
        p = get(ax(2),'pos'); p(4) = p(4) - 0.1; set(ax(2),'pos',p);
        
        axes(ax(3));
        plot(ffrt, ffr, 'o-'); hold on;
        plot(ffrt(isnan(ffr)), min(ffr), 'ro');
        title('Raw (doublec=0.7,freqwt=0.1) ','Fontsize',12);
        
        axes(ax(4));
        plot(ffrt, ffr, 'o-'); hold on;
        plot(ffrt(isnan(ffr)), min(ffr), 'ro');
        title('Raw (doublec=0.7,freqwt=0.1)','Fontsize',12);
        
        for j=1:length(max_F0_jump)
            
             ffr_out = remove_F0_outliers(ffr,sd_tol(k),max_F0_jump(j));
             
             axes(ax(4+j));
             plot(ffrt, ffr_out, 'go-'); hold on;
             plot(ffrt(isnan(ffr_out)), min(ffr_out), 'ro');
             title(['Outliers removed, max-F0-jump = ', num2str(max_F0_jump(j))],'Fontsize',12);
             
        end
        
        linkaxes(ax,'x');
        linkaxes(ax(3:length(ax)),'y');
        
        % Figure title, centered
        set(gcf,'NextPlot','add'); axes;
        h2 = title(['SOME  ', upper(context), ' ', int2str(ind), ' raw  and outlier removed (sd = ', num2str(sd_tol(k)),')'], 'Fontsize', 18);
        set(gca,'Visible','off');
        set(h2,'Visible','on');
        
         % sound(wav,fs);
         
         formatOut = 'mm.dd HHMM';
         timestmp = datestr(now,formatOut);
         filename = [timestmp, ' - ', int2str(ind), ' raw & outlier removed, sd_tol = ', num2str(sd_tol(k)), ', max-F0-jump = ', ...
             regexprep(num2str(max_F0_jump),' *',','),'.eps'];
         print(h, '-depsc', filename);
         
         % pause;
         close(h);
                
    end
    end
%% remove outliers
function [y] = remove_F0_outliers(x,sd_tol,max_F0_jump)
% Bruce's GUI has sd_tol as an input argument; probably some sort of
% standard deviation, but not sure of what. Elsewhere this is hard-coded as
% "2.2" in the first outlier condition.
x = x';

while 1
    
    good_ix = ~isnan(x);
    if sum(good_ix)==0, break; end %contour has no values!
    z = nan(size(x));
    z(good_ix) = zscore(x(good_ix));
    dxp = diff([x(1) x]);  %difference from preceding
    dxn = diff([x x(end)]); %difference from next
    precnan = [false isnan(x(1:end-1))];
    follnan = [isnan(x(2:end)) false];
    
    %outlier conditions
    out_ix1 = abs(z)>sd_tol;                                     %zscore
    out_ix2 = abs(dxn)>max_F0_jump & precnan;                 %big jump to next and preceded by nan
    out_ix3 = abs(dxp)>max_F0_jump & follnan;                 %big jump to preceding and followed by nan
    out_ix4 = (abs(dxn)>max_F0_jump | abs(dxp)>max_F0_jump);  %big jump to either preceding or next
    out_ix5 = ~isnan(x) & precnan & follnan;                  %isolated frame
    
    %exclude outliers
    if sum(out_ix1),        x(out_ix1) = nan;    continue; end
    if sum(out_ix2),        x(out_ix2) = nan;    continue; end
    if sum(out_ix3),        x(out_ix3) = nan;    continue; end
    if sum(out_ix4),        x(out_ix4) = nan;    continue; end
    if sum(out_ix5),        x(out_ix5) = nan;    continue; end
    
    %break if no outliers detected
    if ~sum([out_ix1 out_ix2 out_ix3 out_ix4 out_ix5]), break; end
    
end
y = x;
end


end
