function results=com_ieee8023_93a(config_file, num_fext, num_next, varargin)
% Implementation of Annex 93A of IEEE 802.3
%
% Usage:
% result=com_ieee8023_93a(config_file, num_fext, num_next [, <s4p files>])
% - config_file: XLS file which contains configuration settings (samples in
%   http://www.ieee802.ttorg/3/bj/public/tools/ran_3bj_com_d2p2_02_0813.zip)
% - num_fext: number of FEXT s4p files in the listF
% - num_next: number of NEXT s4p files in the list
% - <s4p_files>: (1+num_fext+num_nefxt) file names. If not supplied, program
%   will ask for each of the files interactively.
%
% This program is intended for the development of channel specifications
% and reflects the work of IEEE P802.3bj, Annex 93A.\
% original proposal may be found at 
% http://www.ieee802.org/3/bj/public/jul12/mellitz_01_0712.pdf in July 2012
% from the following authors and affiliations in 2012.
%      Richard Mellitz, Intel Corporation
%      Charles Moore, Avago Technologies
%      Mike Dudek, QLogic Corporation
%      Mike Peng Li, Altera Corporation
%      Adee Ran, Intel Corporation
%
% This is not an official IEEE document.
%
% Authors and Contributors:
%   Adee Ran 
%   Richard Mellitz 
%   Yasuo Hidaka
% structures:
% chdata(i)
%           i= 1 --> THRU index
%           i= 2, num_fext+1 --> FEXT channel index
%           i= num_fext+2, num_next+num_fext+1
%       base: name of THRU file
%          A: amplitude
%       type: 'THRU', 'NEXT', or 'FEXT'
%        ftr: Rise time frequency
%      fmaxi: max number of frequency points
%      faxis: frequency array [Hz]
%      sdd21: (Htot) rewritten as product of ,vtf based on pkg RL ,tx filter, Rx filter
%      sdd22: differential RL
%      sdd11: differential RL
%     sdd21p: vtf based on pkg RL , this an interim parameter and set to sdd21
%     sdd21f: raw differential IL not filtered use for FD plots
% added output_args.peak_uneq_pulse_mV
% added output_args.cable_loss when "Include PCB" is not 0 in the config file
% added: tap c(-2) c(2) and c(3)
% added: g_DC_HP and f_HP_PZ
% added: new value for "Include PCB" = 2 for cable Rx compliance test only Rx host added 
% added: BREAD_CRUMBS is 1 then a mat file with the structures params and OP is created in the results directory 
% added: T_r_filter_type for RITT testing when IDEAL_TX_TERM is 1: 
% added T_r_meas_point for RITT, if 0, measurement was at tp0, if 1 measurement was tp0a 
% 0 is for is for Gaussian filter and 1 is for a 4th order Bessel-Thompson filter
% fixed INCLUDE_CTLE=0 to really remove from computation
% r161a fixed matlab version issue when for OP.INCLUDE_CTLE=0
% r162 adjusting RITT rise time to Mike Dudek's recommendations also always enable risetime filter if T_r_filter_type=1
% r162 tx and rx package impedance {Zc)
% r162a Gaussian equation corrected 
% r163 cast snr_tx with package test case
% r164 fix pdf for very low noise and lo pass filter enhancements
% r164 add zero gain at nqyist CTLE as in CL12e 
% r165 add simpler congfig command called FORCE_TR (force risetime)
% r200 cm3 and cm4 added cp3 removed
% r200 fixed problem in s21_to_impulse_DC when s parameter have a DC entry
% r200 ILD_FOM updated to EQ93A-55  ERL added 
% r200 improved phase iterpolation for return loss time conversion  
% r200a db = @(x) 20*log10(abs(x)) added so sig processing toolbox won't be required
% r200b Fixed error in bifurcation of Tx/Rx Rd% 
% r200c missed on fixt for interpolation
% r210 new ERL with time gating function
% r211 expertmntal chaael buildibg
%% file_setup
% acquire parsing command string and set up OP control structure. Then read in files
close all; close(findall(0, 'tag', 'TMWWaitbar', '-or', 'tag', 'COM'));
display('This is NOT an official IEEE document.')
display('This is a computation example for exploring COM and ERL for projects like IEEE P802.3 50 Gb/s, 100 Gb/s, and 200 Gb/s Ethernet Task Forcee  P802.3cd ($Revision: 2.11 $)')
display(' with some exploratory extensions and is not normative or official')

%% for review
% close(findall(0, 'Name', 'testing'))
% htest = msgbox('For testing 2.21','testing'); set(htest,'Color',[1 1 .4]);movegui(htest,'northeast');
%%
% need to see what happens for version 8
if verLessThan('matlab', '7.4.1')
    error('Matlab version 7.4 or higher required')
end
% get the first 3 arguments and allow for interactive input.
if ~exist('config_file','var')
    config_file=input('Enter config XLS file or return will just pop a window to ask for the XLS file]:  ','s');
    if isempty(config_file)
        [config_file, config_file_path] = uigetfile({'*.xls'},'INPUT CONFIG FILE .xls');
    end
    if config_file==0
        % cancel - exit gracefully
        return;
    end
    config_file = fullfile(config_file_path, config_file);
end
[param, OP] = read_ParamConfigFile(config_file);
if ~exist('num_fext','var')
    if OP.RX_CALIBRATION
        num_fext=1;
        display('First prompt is for the measured test thru channel and following prompt is for Rx noise path channel')
    else
        if param.tfx == -1 && OP.ERL_ONLY && OP.ERL == 2; % OP.ERL=2 is for package ERL, param.tfx = 1 means fixture time not set and need to be determinind for test fixture channel
            num_fext=1;
            display('First prompt is for the s2p measured data following prompt is for s4p of of the test fixtrure channel')
        elseif ~OP.ERL
            num_fext=input('How many FEXT channels are to be entered? [return means no FEXT] ');
        else
            num_fext=0;
        end
    end
    if isempty(num_fext)==1, num_fext=0; end
end
if ~exist('num_next','var')
    if param.tfx == -1 && OP.ERL_ONLY && OP.ERL == 2; % OP.ERL=2 is for package ERL, param.tfx = 1 means fixture time not set and need to be determinind for test fixture channel
        num_next=0;
    else
        if OP.RX_CALIBRATION
            num_next=0;
        elseif ~OP.ERL_ONLY 
            num_next=input('How many NEXT channels are to be entered? [return means no NEXT] ');
        else
            num_next=0;
        end
    end
    if isempty(num_next)==1, num_next=0; end
end
% Allow string inputs for running compiled version from OS command-line
if ischar(num_fext), num_fext=str2double(num_fext); end
if ischar(num_next), num_next=str2double(num_next); end
xtk=num_fext+num_next; % total number of crosstalk aggressors
param.num_next=num_next;
param.num_fext=num_fext;
param.num_s4p_files=num_fext+num_next+1;
% checking for data when running for rx compliance BBN calibration
if OP.RX_CALIBRATION == 1
    if num_fext ~=1
        h = msgbox('One and only noise path channel is required'); set(h,'Color',[1 .85 0]);
        movegui(h,'northwest')
        if OP.DEBUG ~= 1
            return
        end
    end
    h = msgbox('Please make sure the measured "sigma_RJ", A_DD, and SNR_TX" fields in the config xls file have been modified from the Tx measurement. '); set(h,'Color',[0 1 1]);
    movegui(h,'southeast')
end

if   param.tfx == -1 && OP.ERL_ONLY && OP.ERL == 2;
    if num_fext ~=1
        h = msgbox('One and only test channel is required'); set(h,'Color',[1 .85 0]);
        movegui(h,'northwest')
        if OP.DEBUG ~= 1
            return
        end
    end
    h = msgbox('The test fixture file is use to gate measurements '); set(h,'Color',[0 1 1]);
    movegui(h,'southeast')
end


% create result directory if needed
if ~exist(OP.RESULT_DIR,'dir'); mkdir(OP.RESULT_DIR); end
% allow finite impulse response input rather that s-parameters with
% OP.EXTERNAL = true. However the use_external_IR function is not provided
if ~isempty(varargin) % process case where file names are passed in function call
    if strfind(upper(char(varargin(1))),'EXTERNAL_IR') ~= 0
        OP.EXTERNAL = true;
        OP.GET_FD = 0;
        ir1a= varargin(2);
        ex_var = varargin(3);
        [chdata OP param ]  = use_external_IR(param, OP ,num_fext,num_next,0,ir1a,ex_var);
    else
        OP.EXTERNAL = false;
        if length(varargin) < xtk +1 % check that number of varargin arguments passed is at least number of crosstalk files+1 (thru)
            error('files must include next + fext + a thru');
        end
        %% eveluate any extra arguments as possible modifications of parameters
        extra_args = varargin(xtk+2:end);
        for k=1:2:floor(length(extra_args)/2)*2
            try
                orig_value_is_str = 1;
                orig_value=eval(extra_args{k});
                if ~ischar(orig_value)
                    orig_value_is_str = 0;
                    orig_value=mat2str(orig_value);
                end
            catch eval_err
                if isequal(eval_err.identifier, 'MATLAB:nonExistentField')
                    % trying to modify a nonexistent parameter - probably a
                    % typo. save the user from his error.
                    error('COM:BadExtraParameter', 'Attempted override of a non-existing parameter %s.', extra_args{k});
                else
                    % unexpected condition
                    rethrow(eval_err);
                end
            end
            try
                if orig_value_is_str
                    mod_string = sprintf('%s = ''%s'';', extra_args{k}, extra_args{k+1});
                else
                    mod_string = sprintf('%s = %s;', extra_args{k}, extra_args{k+1});
                end
                eval(mod_string);
                fname=['mod_str' num2str(k)];
                output_args.(fname)=mod_string;
                fprintf('Applied parameter modification: %s (override %s)\n', mod_string, orig_value);
            catch eval_err
                error(eval_err.identifier, 'Error evaluating "%s".', mod_string);
            end
        end
    end
end
%% Parameters computationally defined by values from the settings files
param.ui=1/param.fb;
param.sample_dt = param.ui/param.samples_per_ui;
param.sigma_X=sqrt( (param.levels^2-1)/ (3*(param.levels-1)^2) );
%% size adjust vector parameters which may be entered as one element
if numel(param.C_pkg_board)==1, param.C_pkg_board = param.C_pkg_board*[1 1]; end
if numel(param.C_diepad)==1, param.C_diepad = param.C_diepad*[1 1]; end
if numel(param.R_diepad)==1, param.R_diepad = param.R_diepad*[1 1]; end
if numel(param.CTLE_fp1)==1, param.CTLE_fp1 = param.CTLE_fp1*ones(size(param.ctle_gdc_values)); end
if numel(param.CTLE_fp2)==1, param.CTLE_fp2 = param.CTLE_fp2*ones(size(param.ctle_gdc_values)); end
if numel(param.CTLE_fz)==1, param.CTLE_fz = param.CTLE_fz*ones(size(param.ctle_gdc_values)); end
if numel(param.f_HP)==1, param.f_HP = param.f_HP*ones(size(param.g_DC_HP_values)); end
if numel(param.f_HP_Z)==1, param.f_HP_Z = param.f_HP_Z*ones(size(param.ctle_gdc_values)); end
if numel(param.f_HP_P)==1, param.f_HP_P = param.CTLE_fz*ones(size(param.ctle_gdc_values)); end
% need to adjust to number of cases case 
if numel(param.pkg_Z_c)==1, param.pkg_Z_c = param.pkg_Z_c*[1 1]; end
if numel(param.brd_Z_c)==1, param.brd_Z_c = param.brd_Z_c*[1 1]; end
pkg_sel_vec=ones(1,max(OP.pkg_len_select));
if OP.WC_PORTZ
    if numel(param.a_thru)==1,param.a_thru=param.a_thru*[1 1];end
    if numel(param.a_fext)==1,param.a_fext=param.a_fext*[1 1];end
    if numel(param.a_next)==1,param.a_next=param.a_next*[1 1];end
    if numel(param.SNDR)==1,param.SNDR=param.SNDR*[1 1];end
else
    if numel(param.a_thru)==1,param.a_thru=param.a_thru*pkg_sel_vec;end
    if numel(param.a_fext)==1,param.a_fext=param.a_fext*pkg_sel_vec;end
    if numel(param.a_next)==1,param.a_next=param.a_next*pkg_sel_vec;end
    if numel(param.SNDR)==1,param.SNDR=param.SNDR*pkg_sel_vec;end
end
if ~OP.EXTERNAL
    OP.FIXTURE_CALIBRATION= 0;
    [chdata, param] = get_s4p_files(param, OP, num_fext, num_next, varargin);
end
OP.SAVE_CMD_STR=1;
if OP.SAVE_CMD_STR
 cmd_str = save_cmd_line( config_file, chdata, num_fext,num_next );
 setappdata(0,'cmd_str',cmd_str);
end
%% from here on, multiple package test cases are run. results will be saved separately.
results = cell(size(OP.pkg_len_select));
COM = inf;
min_COM=inf; % reset COM prior to calibration
min_VEO = inf;
sigma_bn=0;
DO_ONCE=true;
low_COM_found = 0;
% at this point only the impulse responses are needed. However vestiges of FD may be intermingled
while (OP.RX_CALIBRATION==1 || DO_ONCE==true) 
    if ~DO_ONCE
        if abs(min_COM - param.pass_threshold)<0.1 || (sigma_bn==0 && min_COM < param.pass_threshold)
            break;
        elseif min_COM > param.pass_threshold
            % increase noise level linearly until low COM found; then perform binary search.
            if low_COM_found
                if OP.sigma_bn_STEP>0 % previous increase too small
                    OP.sigma_bn_STEP = OP.sigma_bn_STEP/2; % gearshift
                else % previously decrease too large
                    OP.sigma_bn_STEP = -OP.sigma_bn_STEP/2; % gearshift and change direction
                end
            end
        else % binary search
            low_COM_found=1;
            if OP.sigma_bn_STEP>0 % previous increase too large
                OP.sigma_bn_STEP = -OP.sigma_bn_STEP/2; % gearshift and change direction
            else % previously decrease too small
                OP.sigma_bn_STEP = OP.sigma_bn_STEP/2; % gearshift
            end
        end
        min_COM = inf; % ignore previous iterations
        min_VEO = inf;
        sigma_bn = sigma_bn + OP.sigma_bn_STEP;
    end
    for package_testcase_i = 1:length(OP.pkg_len_select)
        package_testcase=OP.pkg_len_select(package_testcase_i);
        param.Pkg_len_TX = param.z_p_tx_cases(package_testcase);
        param.Pkg_len_NEXT = param.z_p_next_cases(package_testcase);
        param.Pkg_len_FEXT = param.z_p_fext_cases(package_testcase);
        param.Pkg_len_RX = param.z_p_rx_cases(package_testcase);
        param.package_testcase_i = package_testcase_i;
        if ~OP.EXTERNAL
            %fill in chada with s-parameters
            [chdata, param, SDDch] = read_s4p_files(param, OP, chdata);
            [chdata, param] = process_sxp(param, OP, chdata, SDDch);
            if param.FLAG.S2P, OP.ERL_ONLY =1;end
            if OP.TDR && package_testcase_i == 1
                output_args.Z11=chdata(1).TDR11.avgZport;
                if ~param.FLAG.S2P, output_args.Z22=chdata(1).TDR22.avgZport;end
            end
            if OP.ERL && package_testcase_i == 1
                output_args.ERL11=chdata(1).TDR11.ERL;
                 if ~param.FLAG.S2P,output_args.ERL22=chdata(1).TDR22.ERL;end
                output_args.ERL11RMS=chdata(1).TDR11.ERLRMS;
                 if ~param.FLAG.S2P,output_args.ERL22RMS=chdata(1).TDR22.ERLRMS;end
                
            end
            if OP.ERL_ONLY
                output_args.code_revision='2.21';
                output_args.config_file = config_file;
                output_args.Z_t=param.Z_t;
                fileset_str=str2csv({chdata(1).base});
                output_args.file_names=sprintf('"%s"', fileset_str);
                results = cell(1);
                results{1} = output_args;
                if OP.DISPLAY_WINDOW
                    savefigs(param, OP);
                end
                if OP.CSV_REPORT ==1
                    items = fieldnames(output_args);
                    item_value_strings = cell(size(items));
                    for field_id=1:length(items)
                        field_name=items{field_id};
                        field_value=output_args.(field_name);
                        if ischar(field_value)
                            item_value_strings{field_id}=field_value;
                        elseif numel(field_value)==1
                            item_value_strings{field_id}=num2str(field_value);
                        else
                            item_value_strings{field_id}=sprintf('"%s"', mat2str(field_value));
                        end
                    end
                    
                    header_string = str2csv(items);
                    data_string = str2csv(item_value_strings);
                    if OP.EXTERNAL== false
                        fid = fopen(sprintf('%s%s_case%d_results.csv', OP.RESULT_DIR, chdata(1).base, package_testcase_i),'w');
                        fprintf(fid,'%s\n', header_string);
                        fprintf(fid,'%s\n', data_string);
                        fclose(fid);
                    end
                end  
                break;
            end
        end
        number_of_s4p_files=length(chdata);
        if OP.WC_PORTZ
            A_thru = param.a_thru(param.Tx_rd_sel);
            A_fext = param.a_fext(param.Tx_rd_sel);
            if num_next ~= 0
                A_next = param.a_next(param.Txne_rd_sel);
            end
            param.SNR_TX=param.SNDR(param.Tx_rd_sel);
        else
            A_thru = param.a_thru(package_testcase);
            A_fext = param.a_fext(package_testcase);
            A_next = param.a_next(package_testcase);
            param.SNR_TX=param.SNDR(package_testcase);        
        end
        %% FD processing s-parameter
        % at this point sdd21 responses and faxis (frequency) array are defined
        ICN=0;
        for i=1:number_of_s4p_files
            if isequal(chdata(i).type, 'THRU')
                chdata(i).A=A_thru;
            elseif isequal(chdata(i).type, 'FEXT')
                chdata(i).A=A_fext;
            elseif isequal(chdata(i).type, 'NEXT')
                chdata(i).A=A_next;
            end
        end
        if OP.GET_FD && DO_ONCE
            f2=param.fb;
            MDFEXT_ICN=0; MDNEXT_ICN=0;
            for i=1:number_of_s4p_files
                if i == 2
                    PSXT(1:length(chdata(i).sdd21f))=0;
                    MDFEXT(1:length(chdata(i).sdd21f))=0;
                    MDNEXT(1:length(chdata(i).sdd21f))=0;
                end
                a=find(chdata(i).faxis(:)>f2);
                if isempty(a)
                    f2=chdata(i).faxis(end);
                    index_f2=length(chdata(i).faxis);
                else
                    index_f2=a(1);
                end
                % R is the frequency dependent parameter for the sinc function use in the
                % PWF for ICN
                temp_angle=(param.samples_per_ui*param.sample_dt)*pi.*chdata(i).faxis;
                if(chdata(i).faxis(1)==0)
                    temp_angle(1)=1e-20;% we don't want to divide by zero
                end
                SINC = sin(temp_angle)./temp_angle;
                PWF_data=SINC.^2;
                PWF_trf=(1+(chdata(i).faxis/chdata(i).ftr).^4).^-1;
                %// bw1=2.613126; bw2=3.4142136; bw3=2.613126;
                fr=param.f_r*param.fb;
                PWF_rx=(1+(chdata(i).faxis/fr).^8).^-1;
                PWF_highpass=1;
                %% Equation 93A-57 %%
                PWF=PWF_data.*PWF_trf.*PWF_rx.*PWF_highpass; % power weight function
                % freq delta for integration
                chdata(i).delta_f=chdata(i).faxis(11)-chdata(i).faxis(10);
                % from ba spec, this is basically ICN
                faxis_GHz = chdata(i).faxis/1e9;
                if isequal(chdata(i).type, 'THRU')
                    [ILD chdata(i).fit] = get_ILN(chdata(i).sdd21f(1:index_f2), chdata(i).faxis(1:index_f2));
                    nq_indx=find(chdata(i).faxis >= 1/param.ui/2,1, 'first');
                    chdata(i).fit_ILatNq=-chdata(i).fit(nq_indx);
                    chdata(i).ILatNq=-20*log10(abs(chdata(i).sdd21f(nq_indx)));
                    fit_loss=chdata(i).fit_ILatNq;
                    if OP.include_pcb
                        cable_loss=-20*log10(abs(chdata(i).sdd21_orig(nq_indx)));
                    end
                    total_loss=-20*log10(abs(chdata(i).sdd21p(nq_indx)));
                    Nq_loss=chdata(i).ILatNq;

                    %% Equation 93A-56 %%
                    FOM_ILD=sqrt(chdata(i).delta_f/param.fb*sum( PWF(1:index_f2).*ILD(1:index_f2).^2));
                    if OP.DEBUG
                        if OP.DISPLAY_WINDOW
                            figure(300+package_testcase_i);
                            screen_size=get(0,'ScreenSize');
                            pos = get(gcf, 'OuterPosition');
                            
                            set(gcf, 'Name', [sprintf('%.3gdB IL Channel: ',Nq_loss) 'Raw frequency-domain data'], 'OuterPosition', ...
                                screen_size([3 4 3 4]).*[0 1 0 0] + pos([3 4 3 4]).*[0 -2 1 2] ...
                                - (package_testcase_i-1)*[0 20 0 0]);

                            subplot(3,1,1)
                            title('Losses')
                            plot(faxis_GHz, 20*log10(abs(squeeze(chdata(i).sdd21f))), 'b', 'LineWidth', 3, 'Disp','IL w/o packages')
                            hold on
                            plot(faxis_GHz, 20*log10(abs(squeeze(chdata(i).sdd21))), 'b--', 'Disp','IL with packages')
                            ylim(get(gca, 'ylim'));
                            plot(faxis_GHz(1:index_f2),chdata(i).fit(1:index_f2),'g','Disp','IL fit')
                            plot(faxis_GHz, 20*log10(abs(squeeze(chdata(i).sdd11))),'c','Disp','RL11')
                            plot(faxis_GHz, 20*log10(abs(squeeze(chdata(i).sdd22))),'m','Disp','RL22')
                            subplot(3,1,3)
                            plot(faxis_GHz(1:index_f2), ILD,'Disp','ILD')
                        else
                            display(['Insertion Loss at Nyquist = ', num2str(chdata(i).ILatNq)])
                        end
                    end
                else % NEXT or FEXT
                    if isequal(chdata(i).type, 'FEXT')
                        MDFEXT=sqrt(abs(chdata(i).sdd21f).^2+MDFEXT.^2); % power sum xtk
                        MDFEXT_ICN=sqrt(2*chdata(i).delta_f/param.fb*sum( chdata(i).A^2*PWF(1:index_f2).*abs(MDFEXT(1:index_f2)).^2)); %eq 46
                    elseif isequal(chdata(i).type, 'NEXT')
                        MDNEXT=sqrt(abs(chdata(i).sdd21f).^2+MDNEXT.^2); % power sum xtk
                        MDNEXT_ICN=sqrt(2*chdata(i).delta_f/param.fb*sum( chdata(i).A^2*PWF(1:index_f2).*abs(MDNEXT(1:index_f2)).^2)); %eq 47
                    end
                    PSXT=sqrt((abs(chdata(i).sdd21f)*chdata(i).A).^2+PSXT.^2); % power sum xtk
                    ICN=sqrt(2*chdata(i).delta_f/param.fb*sum( PWF(1:index_f2).*abs(PSXT(1:index_f2)).^2));
                    ICN_test=norm([MDFEXT_ICN MDNEXT_ICN]);
                end
            end  % for loop
            ICN_test=norm([MDFEXT_ICN MDNEXT_ICN]);
            if OP.DEBUG && OP.DISPLAY_WINDOW
                figure(300+package_testcase_i);
                if number_of_s4p_files > 1
                    scale=A_fext/A_thru;% added for crosstalk scaling
                    subplot(3,1,1)
                    hold on
                    plot(faxis_GHz, 20*log10(abs(PSXT*scale)),'r','Disp','PSXTK')
                    icrxi=find(chdata(i).faxis >=param.fb/2,1,'first');
                    subplot(3,1,2)
                    grid on
                    ILtemp=20*log10(abs(chdata(1).sdd21f));
                    IL4ICR=interp1(chdata(1).faxis,ILtemp,chdata(i).faxis);
                    scale=A_fext/A_thru;% added for crosstalk scaling
                    ICR=-20*log10(abs(PSXT*scale))+IL4ICR;
                    semilogx(faxis_GHz, ICR,'Disp', 'ICR')
                    hold on
                    stem(faxis_GHz(icrxi), ICR(icrxi),'g', 'disp', 'f_{Baud}/2')
                end
                subplot(3,1,1)
                title([param.base ' Losses']); ylabel('dB'); xlabel('GHz')
                grid on;  legend show
                subplot(3,1,2)
                title([param.base ' ICR']); ylabel('dB'); xlabel('GHz')
                ylim([0 80])
                xlim([.1 100])
                grid on; %legend show
                subplot(3,1,3)
                title([param.base ' ILD']); ylabel('dB'); xlabel('GHz')
                ylim([-3 3])
                grid on; legend show
            end
        end % get_FD && DO_ONCE
        if DO_ONCE
            % get impulse responses which in interim step between equation for X(f) and
            % H^(k)(t) without TX FFE or CTLE. These will we added later.
            for i=1:number_of_s4p_files
                if OP.INCLUDE_FILTER % apply RX filter
                    %% Equation 93A-20 %%
                    H_r = 1./polyval([1 2.613126 3.414214 2.613126 1], 1i*chdata(i).faxis./(param.f_r*param.fb));
                    chdata(i).sdd21=chdata(i).sdd21.*H_r;
                end
                % if i==1 and else are the same. just do it
                if OP.EXTERNAL == false
                    [chdata(i).uneq_imp_response, ...
                        chdata(i).t, ...
                        chdata(i).causality_correction_dB, ...
                        chdata(i).truncation_dB] = s21_to_impulse_DC(chdata(i).sdd21 ,chdata(i).faxis, param.sample_dt, OP) ;
                    chdata(i).uneq_imp_response=chdata(i).uneq_imp_response*chdata(i).A; % adjust IRx for amplitude
                end
                %------------------------------------------------------------
                % next find Pulse response (SBR) for each channel h^(k)(t)
                if ~OP.DISPLAY_WINDOW && i==1, fprintf('processing COM PDF '); end
                
                chdata(i).uneq_pulse_response=filter(ones(1, param.samples_per_ui), 1, chdata(i).uneq_imp_response);
                
                if OP.DEBUG
                    if OP.DISPLAY_WINDOW && ~OP.RX_CALIBRATION
                        figure(150+package_testcase_i);
                        screen_size=get(0,'ScreenSize');
                        pos = get(gcf, 'OuterPosition');
                        set(gcf, 'OuterPosition', ...
                            screen_size([3 4 3 4]).*[1 1 0 0] + pos([3 4 3 4]).*[-1 -1 1 1] ...
                            - (package_testcase_i-1)*[0 20 0 0]);
                        %movegui(gcf,'northeast')

                        set(gcf, 'Name', sprintf('Case %d PR & PDF - %s', package_testcase_i, chdata(i).base));
                        subplot(2,1,1);  hold on; % all plots on the same axes
                        hp=plot(chdata(i).t, chdata(i).uneq_pulse_response,'Disp', chdata(i).base);
                        % hide thru PR in order to show xtalk in a reasonable
                        % scale. thru is shown in another plot.
                        if isequal(chdata(i).type, 'THRU' ) || OP.RX_CALIBRATION
                            set(hp, 'visible', 'off');
                            set(get(get(hp,'Annotation'),'LegendInformation'), 'IconDisplayStyle','off');
                        end
                        title('Pulse response')
                        ylabel('Volts')
                        xlabel('seconds')
                        
                        recolor_plots(gca);
                    else
                        if param.ndfe~=0
                            fprintf('%s\tUnequalized pulse peak = %.1f mV\n', chdata(i).base, 1000*max(abs(chdata(i).uneq_pulse_response)));
                        end
                    end
                    if OP.ENFORCE_CAUSALITY
                        fprintf('%s\tCausality correction = %.1f dB\n', chdata(i).base, chdata(i).causality_correction_dB);
                    else
                        fprintf('%s\tCausality correction = %.1f dB (not applied)\n', chdata(i).base, chdata(i).causality_correction_dB);
                    end
                    fprintf('%s\tTruncation ratio = %.1f dB\n', chdata(i).base, chdata(i).truncation_dB);
                end
            end
        end %DO_ONCE
        if ~isempty(OP.PHY) % customization hook
            dataout = my_spec(chdata,OP,param);
            try
                items = fieldnames(dataout);
                for field_id=1:length(items)
                    field_name=items{field_id};
                    output_args.(field_name)=dataout.(field_name);
                end
            catch err
                warning(err.identifier, err.message);
            end
        end
        % Determine equalization settings

               
        fom_result = optimize_fom(OP,param, chdata, sigma_bn);
        chdata(1).IRx=fom_result.IR; % IRx has the TX FFE and CTLE applied
        A_s=abs(fom_result.A_s); % this is the "s" in SNR (PAM4 gain in handled in last sections)
        %% Recommended Delta_y no larger than As/1000 or 0.01 mV
        if OP.force_pdf_bin_size
            delta_y = OP.BinSize;
        else
            delta_y = min(A_s/1000, OP.BinSize);
        end
        % the pdf for PAM4 uses the full swing SBR but assigns voltage for the PDF accordingly
        if OP.RX_CALIBRATION, number_of_s4p_files=1; end
        for i=1:number_of_s4p_files
            % get quick PDF  results but only for THRU when in Rx calibration
            if OP.INCLUDE_CTLE==1
                % apply the CTLE that was selected in FOM optimization.
                
%                 ctle_gain = (10^(param.ctle_gdc_values(fom_result.ctle)/20) + 1i*chdata(1).faxis/param.CTLE_fp1(fom_result.ctle)) ./ ...
%                     ((1+1i*chdata(1).faxis/param.CTLE_fp1(fom_result.ctle)).*(1+1i*chdata(1).faxis/param.CTLE_fp2(fom_result.ctle)));
%                 H_low=(10^(param.g_DC_HP_values(fom_result.best_G_high_pass)/20)+  1i*chdata(1).faxis/param.f_HP)./(1 + 1i*chdata(1).faxis/param.f_HP);
%                 H_ctf=H_low.*ctle_gain;
%                 [h_ctf_IR t] = s21_to_impulse_DC(H_ctf ,chdata(1).faxis, 1/param.fb/param.samples_per_ui, OP);
%                 chdata(i).eq_imp_response=filter(h_ctf_IR,1,chdata(i).uneq_imp_response);

                % apply the CTLE that was selected in FOM optimization.
                
                % each  clause gets this but the z adjusted for CL120e
                
                switch param.CTLE_type
                    case 'CL93'
                        chdata(i).eq_imp_response = TD_CTLE(chdata(i).uneq_imp_response, param.fb ...
                            , param.CTLE_fz(fom_result.ctle), param.CTLE_fp1(fom_result.ctle) ...
                            , param.CTLE_fp2(fom_result.ctle), param.ctle_gdc_values(fom_result.ctle) ...
                            , param.samples_per_ui);
                        
                    case 'CL120d'
                        chdata(i).eq_imp_response = TD_CTLE(chdata(i).uneq_imp_response, param.fb ...
                            , param.CTLE_fz(fom_result.ctle), param.CTLE_fp1(fom_result.ctle) ...
                            , param.CTLE_fp2(fom_result.ctle), param.ctle_gdc_values(fom_result.ctle) ...
                            , param.samples_per_ui);
                        
                        chdata(i).eq_imp_response = TD_CTLE(chdata(i).eq_imp_response, param.fb ...
                            , param.f_HP(fom_result.best_G_high_pass), param.f_HP(fom_result.best_G_high_pass) ...
                            , 100e100 , param.g_DC_HP_values(fom_result.best_G_high_pass) ...
                            , param.samples_per_ui);
                    case 'CL120e' % z has been adjusted  for gain
                          chdata(i).eq_imp_response = TD_CTLE(chdata(i).uneq_imp_response, param.fb ...
                            , param.CTLE_fz(fom_result.ctle), param.CTLE_fp1(fom_result.ctle) ...
                            , param.CTLE_fp2(fom_result.ctle), param.ctle_gdc_values(fom_result.ctle) ...
                            , param.samples_per_ui);
                                            
                        chdata(i).eq_imp_response =  TD_CTLE(chdata(i).eq_imp_response,param.fb, ...
                            param.f_HP_Z(fom_result.ctle), param.f_HP_P(fom_result.ctle),1e99, 0, param.samples_per_ui);
                end
            else
                chdata(i).eq_imp_response=chdata(i).uneq_imp_response;
            end

            if isequal(chdata(i).type, 'FEXT') || isequal(chdata(i).type, 'THRU')
                % apply the TXFFE that was selected in FOM optimization.
                upsampled_txffe = zeros(1, param.samples_per_ui*4+1); % start with zeros everywhere
                upsampled_txffe(1+(0:5)*param.samples_per_ui) = fom_result.txffe; % plant the coefficients in the desired locations
                chdata(i).eq_imp_response = filter(upsampled_txffe, 1, chdata(i).eq_imp_response);
            end
            chdata(i).eq_pulse_response=filter(ones(1, param.samples_per_ui), 1, chdata(i).eq_imp_response);
            %chdata(i).result=get_SBR(chdata(i).eq_imp_response, param.samples_per_ui);

            if ~OP.DISPLAY_WINDOW, fprintf('%d ', i); end
            pdf = get_pdf(chdata(i), delta_y, fom_result.t_s, param, OP) ;
            if OP.DEBUG && OP.DISPLAY_WINDOW && ~OP.RX_CALIBRATION
                figure(150+package_testcase_i);
                subplot(2,1,2);
                pdf0.x=pdf.x(pdf.y~=0);
                pdf0.y=pdf.y(pdf.y~=0);
                semilogy(pdf0.x, pdf0.y,'disp', chdata(i).base);
                current_ylim=ylim; ylim([param.specBER/100, current_ylim(2)]);
                hold on; title('PDF')
                recolor_plots(gca);
            end

            chdata(i).pdfr=pdf;
            % reporting
            a=find(cumsum(chdata(i).pdfr.y) >1e-12,1,'first');
            chdata(i).maxquickpdf=(chdata(i).pdfr.y(a));

        end
        if ~OP.DISPLAY_WINDOW, fprintf('\n'); end

        if OP.RX_CALIBRATION
            ctle_gain2 = (10.^(param.ctle_gdc_values(fom_result.ctle)/20) + 1i*chdata(2).faxis/param.CTLE_fz(fom_result.ctle)) ./ ...
                ((1+1i*chdata(2).faxis/param.CTLE_fp1(fom_result.ctle)).*(1+1i*chdata(2).faxis/param.CTLE_fp2(fom_result.ctle)));
            switch param.CTLE_type
                case 'CL93'
                   H_low2=1; 
                case 'CL120d' % this clause uses two gain indexes
                    H_low2=(10.^(param.g_DC_HP_values(fom_result.best_G_high_pass)/20) +  1i*chdata(2).faxis/param.f_HP)./(1 + 1i*chdata(2).faxis/param.f_HP);
                case 'CL120e' % Z1 has been adjusted 
                    H_low2=(1 +  1i*chdata(2).faxis/f_HP_P(fom_result.ctle))./(1 + 1i*chdata(2).faxis/f_HP_Z(fom_result.ctle));
            end
            H_ctf2=H_low2.*ctle_gain2;
            sigma_ne = get_sigma_noise( H_ctf2,  param, chdata, sigma_bn );
        else
            sigma_ne=0;
        end
        sigma_N = fom_result.sigma_N;
        %% Equation 93A-30 %%
        sigma_TX = (param.levels-1)*A_s/param.R_LM*10^(-param.SNR_TX/20);

        %% Equation 93A-41 %%
        sigma_G = norm([param.sigma_RJ*param.sigma_X*norm(fom_result.h_J), sigma_N, sigma_TX]);

        %% Equation 93A-42 %%
        % number of sigmas needed depends on the required BER.
        ber_q = sqrt(2)*erfcinv(2*param.specBER);
        gaussian_noise_pdf = normal_dist(sigma_G, ber_q, delta_y);
        % enable overriding the Q factor of the BBN instrument.
        if OP.force_BBN_Q_factor
            ne_noise_pdf = normal_dist(sigma_ne, OP.BBN_Q_factor, delta_y);
        else
            ne_noise_pdf = normal_dist(sigma_ne, ber_q, delta_y);
        end
        gaussian_noise_pdf = conv_fct(gaussian_noise_pdf, ne_noise_pdf);

        %% p_DD is computed using the procedure defined in 93A.1.7.1 with h(n)=A_DD*h_J(n)
        p_DD = get_pdf_from_sampled_signal(param.A_DD*fom_result.h_J, param.levels, delta_y);

        %% Equation 93A-43
        noise_pdf=conv_fct(gaussian_noise_pdf, p_DD);

        %% Implementation of 93A.1.7.3 combination procedure
        %%  (effectively Equation 93A-44) %%

        % Self-Channel Interference is thru residual result
        sci_pdf = chdata(1).pdfr;
        sci_mxi=find(cumsum(sci_pdf.y)>=param.specBER, 1, 'first');
        thru_peak_interference_at_BER=abs(sci_pdf.x(sci_mxi));
        sci_msi=find(cumsum(sci_pdf.y)>=param.specBER, 1, 'first');
        sci_sigma=abs(sci_pdf.x(sci_msi)/(erfcinv(2*param.specBER)*sqrt(2)));
        if OP.RX_CALIBRATION ==0
            % Co-Channel Interference PDFs (for information only):
            % initialize to deltas
            MDNEXT_cci_pdf = d_cpdf(delta_y, 0, 1);
            MDFEXT_cci_pdf = d_cpdf(delta_y, 0, 1);
            % serially convolve FEXT/NEXT PDFs
            for k=2:number_of_s4p_files
                if isequal(chdata(k).type, 'NEXT')
                    MDNEXT_cci_pdf = conv_fct(MDNEXT_cci_pdf, chdata(k).pdfr);
                else % ... must be FEXT
                    MDFEXT_cci_pdf = conv_fct(MDFEXT_cci_pdf, chdata(k).pdfr);
                end
            end
            
            % find "peaks" of MDNEXT/MDFEXT for reporting
            mdnxi=find(cumsum(MDNEXT_cci_pdf.y)>=param.specBER, 1, 'first');
            MDNEXT_peak_interference=abs(MDNEXT_cci_pdf.x(mdnxi));
            mdfxi=find(cumsum(MDFEXT_cci_pdf.y)>=param.specBER, 1, 'first');
            MDFEXT_peak_interference=abs(MDFEXT_cci_pdf.x(mdfxi));
            
            % Combined crosstalk effect
            cci_pdf = conv_fct(MDFEXT_cci_pdf, MDNEXT_cci_pdf);
            cci_mxi=find(cumsum(cci_pdf.y)>=param.specBER, 1, 'first');
            cci_msi=find(cumsum(cci_pdf.y)>=param.specBER, 1, 'first');
            cci_sigma=abs(cci_pdf.x(cci_msi)/(erfcinv(2*param.specBER)*sqrt(2)));
            crosstalk_peak_interference_at_BER=abs(cci_pdf.x(cci_mxi));
            % combine cci and sci
            isi_and_xtalk_pdf = conv_fct(sci_pdf, cci_pdf);
        else
            % for calibration there is no cci
            isi_and_xtalk_pdf=sci_pdf;
        end
        
        mxi=find(cumsum(isi_and_xtalk_pdf.y)>=param.specBER, 1, 'first');
        peak_interference_at_BER=abs(isi_and_xtalk_pdf.x(mxi));
        
        
        %% Equation 93A-45
        combined_interference_and_noise_pdf = conv_fct(isi_and_xtalk_pdf, noise_pdf);
        
        %% Equation 93A-37
        combined_interference_and_noise_cdf=cumsum(combined_interference_and_noise_pdf.y);
        
        %% The noise and interference amplitude, A_ni, is the magnitude of the value of y0
        %% that satisfies the relationship P(y0) = DER_0
        A_ni_ix=find(combined_interference_and_noise_cdf>param.specBER, 1, 'first');
        A_ni = abs(combined_interference_and_noise_pdf.x(A_ni_ix));
        
        %% Equation 93A-1
        COM=20*log10(A_s/A_ni);
        min_COM = min(min_COM, COM);
        %% Calculation of error propagation and burst probability
        if OP.nburst>0
            % an error burst of length N will cause each of the first N taps tap to mis-correct and create a PAM (2 or
            % 4) noise term - depending on the N'th previous symbol, with double the tap voltage. From this we calculate
            % the probability of staying in error state, i.e. burst of length N+1.
            
            % initialize loop with uncorrelated noise and BER
            error_propagation_noise_pdf{1}=combined_interference_and_noise_pdf; %#ok<AGROW> % PDF for burst of length 1 is the uncorrelated PDF
            
            % Assume an error will occur if the noise excceds the available signal
            % reduced by some dB. reduction is by COM threshold minus Error
            % propagation margin (a positive EP margin reduces uncorrelated error probability
            % below target BER).
            error_threshold =  A_s./10^((param.pass_threshold-OP.COM_EP_margin)/20);
            % Find the probability of this event by integration of the PDF. Use 1e-20 as a floor probabilty if noise PDF isn't wide enough.
            x_error_propagation = find(error_propagation_noise_pdf{1}.x >= error_threshold, 1, 'first');
            if isempty(x_error_propagation)
                p_error_propagation(1) = 1e-20;
            else
                p_error_propagation(1) = sum(error_propagation_noise_pdf{1}.y(x_error_propagation:end));  % uncorrelated BER
            end
            
            sorted_abs_dfe_taps = sort(abs(result.DFE_taps), 'descend');
            for k=2:min(param.ndfe, OP.nburst)
                % (arrays kept to allow tracking during development, though not really needed)
                if OP.use_simple_EP_model
                    post_error_dfe_noise_pdf{k} = get_pdf_from_sampled_signal( 2*A_s*max(sorted_abs_dfe_taps), param.levels, delta_y ); %#ok<AGROW>
                    error_propagation_noise_pdf{k} = conv_fct(error_propagation_noise_pdf{1}, post_error_dfe_noise_pdf{k}); %#ok<AGROW>
                else
                    post_error_dfe_noise_pdf{k} = get_pdf_from_sampled_signal( 2*A_s*sorted_abs_dfe_taps(k-1), param.levels, delta_y ); %#ok<AGROW>
                    error_propagation_noise_pdf{k} = conv_fct(error_propagation_noise_pdf{k-1}, post_error_dfe_noise_pdf{k}); %#ok<AGROW>
                end
                
                % Assume an error will propagate if this noise exceeds the threshold defined above
                x_error_propagation = find(error_propagation_noise_pdf{k}.x >= error_threshold, 1, 'first');
                if isempty(x_error_propagation)
                    p_error_propagation(k) = 1e-20; %#ok<AGROW>
                else
                    p_error_propagation(k) = sum(error_propagation_noise_pdf{k}.y(x_error_propagation:end)); %#ok<AGROW>
                end
            end
            
            % Assume an uncorrelated error will occur if the original noise exceeds
            % the available signal reduced by pass_threhsold dB. Find the probability
            % of this event by partial sum of the PDF.
            % p_uncorrelated_error_i = find(combined_interference_and_noise_pdf.x >= A_s./10^(param.pass_threshold/20), 1, 'first');
            % p_uncorrelated_error = sum(combined_interference_and_noise_pdf.y(p_uncorrelated_error_i:end));
            
            % probability of bursts of different lengths
            p_burst = cumprod(p_error_propagation);
        end
        %% reporting
        output_args.code_revision='2.21';
        %         if exist('cmd_str','var')
        %             output_args.cmd_line=cmd_str;
        %         end
        output_args.config_file = config_file;
        fileset_str=str2csv({chdata.base});
        output_args.file_names=sprintf('"%s"', fileset_str);
        for pkg_params = {'levels', 'Pkg_len_TX', 'Pkg_len_NEXT', 'Pkg_len_FEXT', 'Pkg_len_RX','R_diepad','pkg_Z_c'}
            output_args.(pkg_params{:})= param.(pkg_params{:});
        end
        output_args.baud_rate_GHz=param.fb/1e9;
        output_args.f_Nyquist_GHz = param.fb/2e9;
        output_args.BER=param.specBER;
        output_args.FOM = fom_result.FOM;
        output_args.DFE4_RSS=norm(fom_result.DFE_taps(4:end));
        output_args.DFE2_RSS=norm(fom_result.DFE_taps(2:end));
        output_args.channel_operating_margin_dB=COM;
        output_args.available_signal_after_eq_mV=1000*A_s;
        output_args.peak_uneq_pulse_mV=1000*max(abs(chdata(1).uneq_pulse_response));
        output_args.steady_state_voltage_mV = 1000*sum(fom_result.IR);
        
        if OP.GET_FD
            output_args.FOM_ILD=FOM_ILD;
        end
        if OP.RX_CALIBRATION== 1, output_args.sigma_bn=sigma_bn; end
        output_args.Peak_ISI_XTK_and_Noise_interference_at_BER_mV=1000*A_ni;
        output_args.peak_ISI_XTK_interference_at_BER_mV=1000*peak_interference_at_BER;
        output_args.peak_ISI_interference_at_BER_mV=1000*thru_peak_interference_at_BER;
        output_args.equivalent_ICI_sigma_assuming_PDF_is_Gaussian_mV=sci_sigma*1000;
        
        if OP.RX_CALIBRATION == 0
            output_args.peak_MDXTK_interference_at_BER_mV=1000*crosstalk_peak_interference_at_BER;
            output_args.peak_MDNEXT_interference_at_BER_mV=1000*MDNEXT_peak_interference;
            output_args.peak_MDFEXT_interference_at_BER_mV=1000*MDFEXT_peak_interference;
        end
        output_args.ICN_mV=ICN*1000;
        %             output_args.ICN_test_mV=ICN_test*1000;
        if xtk>0 && OP.RX_CALIBRATION ==0
            output_args.MDNEXT_ICN_92_46_mV=MDNEXT_ICN*1000;
            output_args.MDFEXT_ICN_92_47_mV=MDFEXT_ICN*1000;
            output_args.equivalent_ICN_assuming_Gaussian_PDF_mV=cci_sigma*1000;
        else
            output_args.MDNEXT_ICN_92_46_mV=0;
            output_args.MDFEXT_ICN_92_47_mV=0;
            output_args.equivalent_ICN_assuming_PDF_is_Gaussian_mV=0;
        end
        %output_args.SNR_ISI_XTK_normalized_1_sigma=20*log10(A_s/(peak_interference_at_BER/qfuncinv(param.specBER))); modified by Yasuo Hidaka, 8/7/17
        output_args.SNR_ISI_XTK_normalized_1_sigma=20*log10(A_s/(peak_interference_at_BER/sqrt(2)/erfcinv(2*param.specBER)));
        if OP.GET_FD
            output_args.total_IL_wpkgs_dB_at_Fnq=total_loss;
            output_args.IL_dB_channel_only_at_Fnq=Nq_loss;
            output_args.fitted_IL_dB_at_Fnq = fit_loss;
            if OP.include_pcb
                output_args.cable__assembley_loss=cable_loss;
            end
        end

        if OP.DEBUG
            msg = sprintf('Case %d: z_p=(%d, %d, %d, %d) (TX, RX, NEXT, FEXT):', ...
                package_testcase_i, param.Pkg_len_TX, param.Pkg_len_RX, param.Pkg_len_NEXT, param.Pkg_len_FEXT ...
                );
        else
            msg = sprintf('Case %d:', package_testcase_i );
        end
        VEO_mV = 1000*(A_s-A_ni)*2;
        if OP.DISPLAY_WINDOW && ~OP.RX_CALIBRATION
            % display bathtub curves in one axis per test case.
            if ~OP.COM_CONTRIBUTION_CURVES
                figure_name =  'Voltage bathtub curves';
                fig=findobj('Name', figure_name);
                if isempty(fig), fig=figure('Name', figure_name); end
                figure(fig);
                movegui(fig,'south')
                hax = subplot(length(OP.pkg_len_select), 1, package_testcase_i);
                plot_bathtub_curves( hax ...
                    , A_s ...
                    , sci_pdf ...
                    , cci_pdf ...
                    , isi_and_xtalk_pdf ...
                    , noise_pdf ...
                    , combined_interference_and_noise_pdf ...
                    , delta_y ...
                    );
                set(hax, 'tag', 'BTC');
                title(hax, sprintf('case %d VBC: %s ', package_testcase_i, regexprep([chdata(1).base,' '],'_',' ')));
                ylim(hax, [param.specBER/10 1]);
                % show BER target line
                hp=plot(get(hax, 'xlim'), param.specBER*[1 1], 'r:');
                set(get(get(hp,'Annotation'),'LegendInformation'), 'IconDisplayStyle','off');
            else
                figure_name =  'COM Contributions (Rough Allocations)';
                fig=findobj('Name', figure_name);
                if isempty(fig), fig=figure('Name', figure_name); end
                figure(fig);
                movegui(fig,'south')
                hax = subplot(length(OP.pkg_len_select), 1, package_testcase_i);
                
                plot_pie_com( hax ...
                    , A_s ...
                    , sci_pdf ...
                    , cci_pdf ...
                    , isi_and_xtalk_pdf ...
                    , noise_pdf ...
                    , combined_interference_and_noise_pdf ...
                    , delta_y, param...
                    );
                set(hax, 'tag', 'BTC');
                title(hax, sprintf('case %d rough COM impact: %s ', package_testcase_i, regexprep([chdata(1).base,' '],'_',' ')));
            end
            h=findall(0, 'Name', 'COM results');
            if isempty(h)
                msgtext = cell(1, length(OP.pkg_len_select));
                msgcolor = 'g';
            else
                msgtext=get(findobj(h, 'type', 'text'), 'string');
                msgcolor = get(h, 'color');
                close(h); % will be recreated
            end
            switch param.PHY_type
                case 'C2M'
                    if VEO_mV >= param.Min_VEO && VEO_mV <= param.Max_VEO
                        msgtext{package_testcase_i}=sprintf('%s: EH = %.3f mV (pass)\n', ...
                            msg, VEO_mV);
                    else
                        msgtext{package_testcase_i}=sprintf('%s: EH = %.3f mV (FAIL)\n', ...
                            msg, VEO_mV);
                        msgcolor = 'r';
                    end
                case 'C2C'
                    if COM >= param.pass_threshold
                        msgtext{package_testcase_i}=sprintf('%s: COM = %.3f dB (pass)\n', ...
                            msg, COM);
                    else
                        msgtext{package_testcase_i}=sprintf('%s: COM = %.3f dB (FAIL)\n', ...
                            msg, COM);
                        msgcolor = 'r';
                    end
            end
            close(findall(0, 'tag', 'COM'));
            h=msgbox(msgtext, 'COM results');
            set(h, 'color', msgcolor, 'tag', 'COM');
            movegui(h, 'center');
        else % no windows
            display(['max noise at BER = ' num2str(peak_interference_at_BER)])
            display(['signal after eq = ' num2str(A_s/(param.levels-1))])
            switch param.PHY_type
                case 'C2C'
                    if COM >= param.pass_threshold
                        fprintf('%s <strong> PASS ... COM = %.3f dB</strong>\n', msg, COM);
                    else
                        fprintf('%s <strong> FAIL ... COM = %.3f dB</strong>\n',  msg, COM);
                    end
                case 'C2M'
                    if VEO_mV >= param.Min_VEO && VEO_mV <= param.Max_VEO
                        fprintf('%s <strong> PASS ... EH = %.3f mV</strong>\n', msg, VEO_mV); 
                    else                        
                        fprintf('%s <strong> FAIL ... EH = %.3f mV</strong>\n',  msg, VEO_mV);
                    end
            end
        end
        switch param.CTLE_type
            case 'CL93'
                output_args.CTLE_zero_poles=[param.CTLE_fz(fom_result.ctle) param.CTLE_fp2(fom_result.ctle) param.CTLE_fp1(fom_result.ctle)];
                output_args.CTLE_DC_gain_dB=param.ctle_gdc_values(fom_result.ctle);
            case 'CL120d'
                output_args.CTLE_zero_poles=[param.CTLE_fz(fom_result.ctle) param.CTLE_fp2(fom_result.ctle) param.CTLE_fp1(fom_result.ctle)];
                output_args.CTLE_DC_gain_dB=param.ctle_gdc_values(fom_result.ctle);
                output_args.g_DC_HP=param.g_DC_HP_values(fom_result.best_G_high_pass);
                output_args.HP_poles_zero=param.f_HP(fom_result.best_G_high_pass);
            case 'CL120e'
                output_args.CTLE_zero_poles=[param.CTLE_fz(fom_result.ctle) param.f_HP_Z(fom_result.ctle) param.CTLE_fp2(fom_result.ctle) param.CTLE_fp1(fom_result.ctle) param.f_HP_P(fom_result.ctle)];
                output_args.CTLE_DC_gain_dB=param.ctle_gdc_values(fom_result.ctle);
        end
        output_args.TXLE_taps=fom_result.txffe;
        output_args.DFE_taps=fom_result.DFE_taps;
        if OP.nburst>0
            output_args.error_propagation_probability = p_error_propagation;
            output_args.burst_probabilities = p_burst;
        end
        %% making mat file
        if ~OP.EXTERNAL
            if(OP.DEBUG)
                save (sprintf('%s%s_case%d_results.mat', OP.RESULT_DIR, chdata(1).base, package_testcase_i), ...
                    'output_args');
            end
        else %??
            output_args.baud_rate_GHz=param.fb/1e9;
            output_args.sigma_rj=param.sigma_RJ;
            output_args.Add_noise=param.A_DD;
        end
        min_VEO = min(min_VEO, VEO_mV);
        output_args.VEO_normalized = (A_s-A_ni)/A_s;
        output_args.VEC_dB = -20*log10(output_args.VEO_normalized);
        output_args.VEO_mV = 1000*(A_s-A_ni)*2;
        output_args.COM_dB=COM;
        %% making csv file
        if OP.CSV_REPORT ==1
            items = fieldnames(output_args);
            item_value_strings = cell(size(items));
            for field_id=1:length(items)
                field_name=items{field_id};
                field_value=output_args.(field_name);
                if ischar(field_value)
                    item_value_strings{field_id}=field_value;
                elseif numel(field_value)==1
                    item_value_strings{field_id}=num2str(field_value);
                else
                    item_value_strings{field_id}=sprintf('"%s"', mat2str(field_value));
                end
            end
            
            header_string = str2csv(items);
            data_string = str2csv(item_value_strings);
            if OP.EXTERNAL== false
                fid = fopen(sprintf('%s%s_case%d_results.csv', OP.RESULT_DIR, chdata(1).base, package_testcase_i),'w');
                fprintf(fid,'%s\n', header_string);
                fprintf(fid,'%s\n', data_string);
                fclose(fid);
            end
        end
        results{package_testcase_i} = output_args;
        if nargout==0
            fprintf('<strong>--- Testcase %d results ---</strong>\n', package_testcase_i);
            disp(output_args)
        end
        % package test cases
        
        if OP.DEBUG && OP.DISPLAY_WINDOW && OP.RX_CALIBRATION==0
            btc_axes = findobj('tag', 'BTC');
            if ~isempty(btc_axes), linkaxes(btc_axes, 'x'); end
            eqe_axes = findobj('tag', 'EQE');
            if ~isempty(eqe_axes), linkaxes(eqe_axes, 'xy'); end
        end
    end
    if ~OP.DEBUG
        % display bottom line results, minimum COM across test cases
        if OP.DISPLAY_WINDOW
            close(findall(0, 'tag', 'COM'));
            switch param.PHY_type
                case 'C2C'
                    if min_COM >= param.pass_threshold
                        h = msgbox(sprintf('PASS ... COM = %.3f dB', min_COM));
                        set(h,'Color','g', 'tag', 'COM');
                    else
                        h = msgbox(sprintf('FAIL ... COM = %.3f dB', min_COM));
                        set(h,'Color','r', 'tag', 'COM');
                    end
                case 'C2M'
                    if min_VEO >= param.VEO
                        h = msgbox(sprintf('PASS ... VEO = %.3f mV', min_VEO));
                        set(h,'Color','g', 'tag', 'COM');
                    else
                        h = msgbox(sprintf('FAIL ... VEO = %.3f mV', min_VEO));
                        set(h,'Color','r', 'tag', 'COM');
                    end
                case 'C2M'
            end
        end
    end
    if ~OP.ERL_ONLY
        switch param.PHY_type
            case 'C2C'
                if min_COM >= param.pass_threshold
                    fprintf('<strong> PASS ... COM = %.3f dB</strong>\n', min_COM);
                else
                    fprintf('FAIL ... COM = %.3f dB\n', min_COM);
                end
            case 'C2M'
                if VEO_mV >= param.Min_VEO && VEO_mV <= param.Max_VEO
                    fprintf('<strong> PASS ... EH = %.3f mV</strong>\n', min_VEO);
                else
                    fprintf('FAIL ... EH = %.3f mV\n', min_VEO);
                end
        end
        if OP.ERL
            fprintf('<strong> ERL11 = %.3f dB:  ERL22 = %.3f dB:</strong>\n', output_args.ERL11, output_args.ERL22);
            fprintf('<strong> ERL11RMS = %.3f dB:  ERL22RMS = %.3f dB:</strong>\n', output_args.ERL11RMS, output_args.ERL22RMS);
        end
    else
        if ~param.FLAG.S2P
            fprintf('<strong> ERL11 = %.3f dB:  ERL22 = %.3f dB:</strong>\n', output_args.ERL11, output_args.ERL22);
            fprintf('<strong> ERL11RMS = %.3f dB:  ERL22RMS = %.3f dB:</strong>\n', output_args.ERL11RMS, output_args.ERL22RMS);
        else
            fprintf('<strong> ERL11 = %.3f dB:</strong>\n', output_args.ERL11);
            fprintf('<strong> ERL11RMS = %.3f dB</strong>\n', output_args.ERL11RMS);

        end
    end
    %%
    if OP.RX_CALIBRATION ==1
        display ([' LOOP with sigma_bn = ' num2str(sigma_bn) ' performed with COM = ' num2str(min_COM) ])
    end
    DO_ONCE=false;
end

if OP.DISPLAY_WINDOW
    savefigs(param, OP);
end

if OP.RX_CALIBRATION==1
    fprintf ('Set Rx calibration noise rms voltage to %g mV\n', sigma_bn*1000);
    if OP.DISPLAY_WINDOW
        message=sprintf('Set Rx calibration noise rms voltage to %g mV.',sigma_bn*1000);
        hlast = msgbox(message,'sigma_bn','help');
        set(hlast,'Color','y', 'tag', 'COM');
    end
end
if OP.BREAD_CRUMBS
    save( [OP.RESULT_DIR 'param_n_OP_' param.base '.mat'],'param','OP')
end

if length(results)==1, results = results{1}; end

if isdeployed
    if OP.exit_if_deployed
        quit
    end
end
%%
%--------------------------------------------------------------------------
%--------------- Helper functions -----------------------------------------
%--------------------------------------------------------------------------
function [voltage, t_base, causality_correction_dB, truncation_dB] = ...
    s21_to_impulse_DC(IL, freq_array, time_step, OP)
% Creates a time-domain impulse response from frequency-domain IL data.
% IL does not need to have DC but a corresponding frequency array
% (freq_array) is required.
%
% Causality is imposed using the Alternating Projections Method. See also:
% Quatieri and Oppenheim, "Iterative Techniques for Minimum Phase Signal
% Reconstruction from Phase or Magnitude", IEEE Trans. ASSP-29, December
% 1981 (http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1163714)

ILin=IL;
fmax=1/time_step/2;
freq_step=(freq_array(3)-freq_array(2))/1;
fout=0:1/round(fmax/freq_step)*fmax:fmax;
IL=interp_Sparam(ILin,freq_array,fout, ...
    OP.interp_sparam_mag, OP.interp_sparam_phase,OP);
IL_nan = find(isnan(IL));
for in=IL_nan
    IL(in)=IL(in-1);
end

IL = IL(:);
% add padding for time steps
IL_symmetric = [IL(1:end-1);0; flipud(conj(IL(2:end-1)))];
impulse_response = real(ifft(IL_symmetric));
L = length(impulse_response);
t_base = (0:L-1)/(freq_step*L);

original_impulse_response=impulse_response;
% Correct non-causal effects frequently caused by extrapolation of IL
% Assumption: peak of impulse_response is in the first half, i.e. not anti-causal
abs_ir=abs(impulse_response);
a = find(abs_ir(1:L/2) > max(abs_ir(1:L/2))*OP.EC_PULSE_TOL);
start_ind = a(1);

err=inf;
while ~all(impulse_response==0)
    impulse_response(1:start_ind)=0;
    impulse_response(floor(L/2):end)=0;
    IL_modified=abs(IL_symmetric).*exp(1j*angle(fft(impulse_response)));
    ir_modified = real(ifft(IL_modified));
    delta = abs(impulse_response-ir_modified);
    
    err_prev = err;
    err=max(delta)/max(impulse_response);
    if err<OP.EC_REL_TOL || abs(err_prev-err)<OP.EC_DIFF_TOL
        break;
    end
    
    impulse_response=ir_modified;
end

causality_correction_dB=20*log10(norm(impulse_response-original_impulse_response)/norm(impulse_response));

if ~OP.ENFORCE_CAUSALITY
    impulse_response = original_impulse_response;
end
% truncate final samples smaller than 1e-3 of the peak
ir_peak = max(abs(impulse_response));
ir_last  = find(abs(impulse_response)>ir_peak*OP.impulse_response_truncation_threshold, 1, 'last');

voltage = impulse_response(1:ir_last);
t_base = t_base(1:ir_last);

truncation_dB=20*log10(norm(impulse_response(ir_last+1:end))/norm(voltage));


function [Sout] = interp_Sparam(Sin,fin,fout, ...
    opt_interp_Sparam_mag, opt_interp_Sparam_phase,OP)
% Sout = interp_Sparam(Sin,fin,fout)
%
% Interpolate S-parameters Sin from frequency grid fin to frequency grid
% fout.

if ( fin(end)<fout(end) )
    %    warning('Channel high frequencies extrapolation might be inaccurate!');
end

H_mag = abs(Sin);
H_mag(H_mag<eps)=eps; % handle ill cases...
H_ph = unwrap(angle(Sin));
% For long delay channels, the result can turn anti-causal if frequency step is too coarse. Don't let the
% user ignore that.
if mean(diff(H_ph))>0
    if OP.DEBUG
        warning('Anti-causal response found. Finer frequency step is required for this channel');
    else
            error('Anti-causal response found. Finer frequency step is required for this channel');
    end
end

%opt_interp_Sparam_mag='linear_trend_to_DC';
switch opt_interp_Sparam_mag
    case 'linear_trend_to_DC'
        if -iscolumn(H_mag), H_mag=H_mag.';end
        if -iscolumn(fin), fin=fin.';end
        fin_x=fin;
        H_mag_x=H_mag(:);
        if fin(1)>0
            p=polyfit(fin(1:10), H_mag(1:10), 1);
            dc_trend_val=polyval(p, 0);
            fin_x=[0, fin_x];
            H_mag_x = [dc_trend_val; H_mag_x];
        end
        if fin(end)<fout(end)
            warn_state=warning('off', 'MATLAB:polyfit:RepeatedPointsOrRescale');
            mid_freq_ind=round(length(fin)/2);
            p=polyfit(fin(mid_freq_ind:end), H_mag(mid_freq_ind:end), 1);
            warning(warn_state);
            hf_trend_val=polyval(p, fout(end));
            if hf_trend_val>H_mag(end)
                hf_trend_val=H_mag(end);
            elseif hf_trend_val<eps
                hf_trend_val=eps;
            end
            fin_x=[fin_x, fout(end)];
            H_mag_x = [H_mag_x; hf_trend_val];
        end
        H_mag_i = interp1(fin_x, H_mag_x, fout, 'linear', 'extrap');
    case 'trend_to_DC'
        % extrapolate to trend value at DC.
        if -iscolumn(H_mag), H_mag=H_mag.';end
        if -iscolumn(fin), fin=fin.';end     
        fin_x=fin;
        H_mag_x=H_mag;
        if fin(1)>0
            p=polyfit(fin(1:10), log10(H_mag(1:10)), 1);
            dc_trend_val=10^polyval(p, 0);
            fin_x=[0, fin_x];
            H_mag_x = [dc_trend_val H_mag_x];
        end
        if fin(end)<fout(end)
            warn_state=warning('off', 'MATLAB:polyfit:RepeatedPointsOrRescale');
            mid_freq_ind=round(length(fin)/2);
            p=polyfit(fin(mid_freq_ind:end), log10(H_mag(mid_freq_ind:end)), 1);
            warning(warn_state);
            hf_trend_val=10^polyval(p, fout(end));
            if hf_trend_val>H_mag(end)
                hf_trend_val=H_mag(end);
            end
            fin_x=[fin_x, fout(end)];
            H_mag_x = [H_mag_x hf_trend_val];
        end
        H_mag_i = 10.^interp1(fin_x,log10(H_mag_x),fout,'linear', 'extrap');
    case 'extrap_to_DC_or_zero'
        % same as extrap_to_DC but detect AC-coupled channels and
        % extrapolate them to 0.
        if fin(1)>0 && 20*log10(H_mag(1))<-20
            % assume AC coupling, with 0 at DC
            H_mag_i = 10.^interp1([0, fin],[-100; log10(H_mag)],fout(fout<=fin(end)),'linear', 'extrap');
        else
            H_mag_i = 10.^interp1(fin, log10(H_mag), fout(fout<=fin(end)),'linear', 'extrap');
        end
        H_mag_i(fout>fin(end)) = H_mag(end);
    case 'extrap_to_DC'
        % first extrapolate down to DC, then use highest available frequency
        % for higher frequencies
        H_mag_i = 10.^interp1(fin,log10(H_mag),fout(fout<=fin(end)),'linear', 'extrap');
        H_mag_i(fout>fin(end)) = H_mag(end);
    case 'old'
        H_mag_i = interp1(fin,H_mag,fout,'linear','extrap');
    otherwise
        error('COM:Extrap:InvalidOption', 'opt_interp_Sparam_mag valid values are "old", "extrap_to_DC"');
end

H_ph_i = interp1(fin,H_ph,fout,'linear', 'extrap');

%opt_interp_Sparam_phase='trend_and_shift_to_DC';
%opt_interp_Sparam_phase='interp_cubic_to_dc_linear_to_inf';
switch opt_interp_Sparam_phase
    case 'old'
        H_ph_i = H_ph_i-H_ph_i(1);
    case 'zero_DC'
        H_ph_i(1) = 0;
    case 'interp_to_DC'
        if fin(1) ~= 0
            H_ph_i = interp1([0; fin(:)], [0; H_ph(:)], fout, 'linear', 'extrap');
        end
    case 'extrap_cubic_to_dc_linear_to_inf'
        if fin(1) ~= 0
            % estimate low frequency group delay
            group_delay = -diff(H_ph(:))./diff(fin(:));
            low_freq_gd = group_delay(1:50);
            %  calculate trend, throwing away outliers
            m = median(low_freq_gd); sigma = std(low_freq_gd);
            lf_trend = mean(low_freq_gd(abs(low_freq_gd-m)<sigma));
            % correct outliers in first 10 phase samples
            for k=10:-1:1
                H_ph(k) = H_ph(k+1) + lf_trend*(fin(k+1)-fin(k));
            end
            H_ph_cubic = interp1(fin, H_ph, fout, 'pchip', 'extrap');
            H_ph_linear = interp1(fin, H_ph, fout, 'linear', 'extrap');
            % modification - trend to inf
            if (1)
                high_freq_gd = group_delay(end-50:end);
                %  calculate trend, throwing away outliers
                m = median(high_freq_gd); sigma = std(high_freq_gd);
                hf_trend = -mean(high_freq_gd(abs(high_freq_gd-m)<sigma));
                hf_extrap_range = find(fout>fin(end));
                last_data_sample = hf_extrap_range(1)-1;
                H_ph_linear(hf_extrap_range) = H_ph_linear(last_data_sample) + (fout(hf_extrap_range)-fout(last_data_sample))*hf_trend;
%                for k=hf_range
%                    H_ph_linear(k) = H_ph_linear(k-1) + hf_trend*(fout(k)-fout(k-1));
%                end
            end

            [UNUSED_OUTPUT, indx] = min(abs(H_ph_cubic-H_ph_linear)); %#ok<ASGLU>
            H_ph_i=H_ph_cubic;
            H_ph_i(indx:end) = H_ph_linear(indx:end);
        end
    case 'interp_and_shift_to_DC'
        if fin(1) ~= 0
            dc_phase_trend = H_ph(1)-(H_ph(2)-H_ph(1))/(fin(2)-fin(1))*fin(1);
            H_ph_i = interp1([0; fin(:)], [0; H_ph(:)-dc_phase_trend], fout, 'linear', 'extrap');
        end
    case 'trend_and_shift_to_DC'
        % estimate low frequency group delay
        group_delay = -diff(H_ph(:))./diff(fin(:));
        low_freq_gd = group_delay(1:50);
        %  calculate trend, throwing away outliers
        m = median(low_freq_gd); sigma = std(low_freq_gd);
        lf_trend = mean(low_freq_gd(abs(low_freq_gd-m)<sigma));
        fin_x=fin;
        H_ph_x=H_ph(:);
        if fin(1) ~= 0
            % correct outliers in first 10 phase samples
            for k=10:-1:1
                H_ph(k) = H_ph(k+1) + lf_trend*(fin(k+1)-fin(k));
            end

            % shift all phase data so that DC extrapolation to 0 follows trend
            dc_phase_trend = H_ph(1)+lf_trend*(fin(1)-0);
            fin_x=[0, fin_x];
            H_ph_x=[0; H_ph(:)-dc_phase_trend];
        end
        % Modification: extrapolate using trend. (interp1 with "extrap" extrapolates using just
        % the last two samples, so noise can create an inverted slope and
        % non-causal response).
        if fout(end)>fin(end)
            group_delay = -diff(H_ph_x(:))./diff(fin_x(:));
            % p=polyfit(fin_x', H_ph_x, 1);
            hf_phase_trend = H_ph_x(end)-median(group_delay)*(max(fout)-max(fin_x));
            % hf_phase_trend=polyval(p,max(fout));
            fin_x=[fin_x, fout(end)];
            H_ph_x=[H_ph_x; hf_phase_trend];
        end
        H_ph_i = interp1(fin_x, H_ph_x, fout, 'linear', 'extrap');

    otherwise
        error('COM:Extrap:InvalidOption', ...
            'debug_interp_Sparam valid values are "old", "zero_DC", "interp_to_DC", "interp_and_shift_to_DC", "trend_and_shift_to_DC", "interp_cubic_to_dc_linear_to_inf"');
end
H_i = H_mag_i.*exp(1j*H_ph_i);
Sout=H_i;

function [data, SDD, SDC] = read_p4_s4params(infile, plot_ini_s_params, plot_dif_s_params, ports,OP)
%% FUNCTION :: read_sp4_sparams
%
% Description
%   Read the fid of single-ended 4-port complex S-parameters
%   in Touchstone format 'file' and convert to the internal
%   format using the port transform 'ports'
%
%   Created by Mike Y. He
%   April 22, 2005
%
%   Reused some code from
%   Anthony Sanders, Alex Deas, Bob Davidov (24 January 2005)
%   for touchstone 4-port S-matrix import.
%
%   Modified (2012-July-27) by Ken Young to match current indexing scheme and
%   optimized for quicker parameter matching and parsing. also, separated out
%   the plotting algorithms into their own sub-function routines
%
% Input Variables (required)
%   file                -- The s4p file to be read and converted
%   plot_ini_sparams    -- Plot the initial s-parameter information. For debugging purposes
%   plot_diff_sparams   -- Plot the differential s-parameter information. For debugging purposes
%   swap                -- Re-order the port layout
%
% Output/Return Variables
%   sch         -- the data matrix contain the network parameter data points
%   schFreq     -- the frequency vector of the network parameter data points
%   sdc         -- the differential in/common-mode out s-parameter data matrix
%   sdd         -- the differential in/differential out s-parameter data matrix
%
% Local Variables (scoped locally to current function only)
%   fid                 -- the opened file being processed
%   fileDataLine        -- the current line in the fid
%   optionLine          -- the options line in the s4p file
%   ports               -- the order of the ports being processed
%   idxNetParams        -- the position (line) of the network parameter data points in the s4p file
%   netParamDataPoints  -- the network parameter data points of the current line of the fid
%   schFreq             -- the frequency points in the s4p file stored as a separate vector
%   sch                 -- the network parameter data points from the s4p file converted in matrix form
%   S                   --
%   T                   --
%   W                   --
%   D                   --
%   sdd                 -- the differential in/differential out form of network parameter data points
%   sdc                 -- the differential in/common-mode out form of network parameter data points
%   plot_ini_sparams    -- Plot the initial s-parameter information. For debugging purposes
%   plot_diff_sparams   -- Plot the differential s-parameter information. For debugging purposes
%

% open the fid
fid = fopen(infile, 'r');

data_file = textscan(fid, '%s', 'delimiter', sprintf('\n'));
data_file = data_file{1,1}; % removes a dimension layer of cell array depth

fclose(fid);

% backwards compatibility settings. can be removed in updated code.
if ~exist('OP', 'var'); OP.DISPLAY_WINDOW = true; end
if isempty(ports); ports = [1 3 2 4]; end % default order normally used.

% adjust ports to maintain the meaning [in1, in2 , out1, out2] when one
% pair is reversed.
ports_adj=ports; for k=1:4, ports_adj(k)=find(ports==k); end; ports=ports_adj;

if OP.DISPLAY_WINDOW
    set(0,'defaulttextinterpreter','none') % prevents subscripting character in displayed messages
    hMsgBox = msgbox(infile, 'Reading S-Parameter File'); % display a progress bar for reading the s-parameter file(s)
end

% parse the file array until the s-parameters options line is found
idxLine = 1;
while true
    if regexp(data_file{idxLine}, '^#', 'once')
        % tests the option fileDataLine cell array to see if it contains...
        if ~isempty(regexpi(data_file{idxLine}, '\sGHz')); freq_scale = 1e9; end
        if ~isempty(regexpi(data_file{idxLine}, '\sMHz')); freq_scale = 1e6; end
        if ~isempty(regexpi(data_file{idxLine}, '\sKHz')); freq_scale = 1e3; end
        if ~isempty(regexpi(data_file{idxLine}, '\sHz'));  freq_scale = 1e0; end
        if ~isempty(regexpi(data_file{idxLine}, 'RI')) && ~isempty(regexpi(data_file{idxLine}, 'R'))
            format = 'complex';
        end
        if ~isempty(regexpi(data_file{idxLine}, 'MA')) && ~isempty(regexpi(data_file{idxLine}, 'R'))
            format = 'linpolar';
        end
        if ~isempty(regexpi(data_file{idxLine}, 'DB')) && ~isempty(regexpi(data_file{idxLine}, 'R'))
            format = 'dBpolar';
        end

        % if the format or freq_scale are not set then exit with error
        if isempty(freq_scale)
            error( 'read_s4p:MissingOptionsParameter', ['\n\t THE FREQUENCY SCALE PARAMETER INFORMATION IS MISSING' ...
                'FROM THE OPTIONS LINE IN THE S4P FILE']);
        elseif isempty(format)
            error( 'read_s4p:MissingOptionsParameter', ['\n\t THE FORMAT PARAMETER INFORMATION IS MISSING' ...
                'FROM THE OPTIONS LINE IN THE S4P FILE']);
        end

        idxLine = idxLine+1; %step one more line before exiting loop
        break
    end
    idxLine = idxLine+1;
end

col_1=1; col_2=2; col_3=3; col_4=4; % This is to maintain current indexing compatibility as well as readability
schFreqAxis = zeros(1,int8((length(data_file)-idxLine)/4)); % preallocated array length for memory and time purposes
idxNetParams = 1; % each 4-port network parameter data point set will receive an index

freqCounter = 0; % variable that will be used to capture number of freq points.

while true
    if isempty(data_file{idxLine}) || ~isempty(regexp(data_file{idxLine}, '^!', 'once'))
        idxLine = idxLine+1;
    else
        for row=1:4
%             dataPointsLine = eval (['[' data_file{idxLine} ']']);
            dataPointsLine= sscanf(data_file{idxLine}, '%f'); % Close V2015 memory leak - Adam Gregory
            if mod(length(dataPointsLine),8)==1
                schFreqAxis(idxNetParams) = dataPointsLine(1) * freq_scale;
                dataPointsLine(1) = []; % removes the frequency value from the current data fileDataLine cell array
                % so that only the network parameter data is in the array
            end
            if mod(length(dataPointsLine),8)~=0  % sanity check to make sure the network parameter data point line structure is correct
                msg = sprintf(['\n'...
                    '\n\t ******************************************************' ...
                    '\n\t ** There is an error in the s4p network parameter   **' ...
                    '\n\t ** data point information structure at line %G.'        ...
                    '\n\t **              ENDING PROGRAM                      **' ...
                    '\n\t ******************************************************' ...
                    '\n\n\n\t Press any key to exit (an error will be displayed). \n\n']);
                disp(msg)
                pause on;
                pause;
                pause off;
                error( 'COM:reads_4p:SparamFileDataError', 'INCORRECT S4P NETWORK PARAMETER DATA POINT STRUCTURE');
            else
                switch format
                    case 'complex'
                        sch(idxNetParams, ports(row), ports(col_1)) = dataPointsLine(1) + 1i*dataPointsLine(2); %#ok<AGROW>
                        sch(idxNetParams, ports(row), ports(col_2)) = dataPointsLine(3) + 1i*dataPointsLine(4); %#ok<AGROW>
                        sch(idxNetParams, ports(row), ports(col_3)) = dataPointsLine(5) + 1i*dataPointsLine(6); %#ok<AGROW>
                        sch(idxNetParams, ports(row), ports(col_4)) = dataPointsLine(7) + 1i*dataPointsLine(8); %#ok<AGROW>
                    case 'linpolar'
                        sch(idxNetParams, ports(row), ports(col_1)) = dataPointsLine(1) * exp(1i*dataPointsLine(2)*(pi/180)); %#ok<AGROW>
                        sch(idxNetParams, ports(row), ports(col_2)) = dataPointsLine(3) * exp(1i*dataPointsLine(4)*(pi/180)); %#ok<AGROW>
                        sch(idxNetParams, ports(row), ports(col_3)) = dataPointsLine(5) * exp(1i*dataPointsLine(6)*(pi/180)); %#ok<AGROW>
                        sch(idxNetParams, ports(row), ports(col_4)) = dataPointsLine(7) * exp(1i*dataPointsLine(8)*(pi/180)); %#ok<AGROW>
                    case 'dBpolar'
                        sch(idxNetParams, ports(row), ports(col_1)) = 10^(dataPointsLine(1)/20) * exp(1i*dataPointsLine(2)*(pi/180)); %#ok<AGROW>
                        sch(idxNetParams, ports(row), ports(col_2)) = 10^(dataPointsLine(3)/20) * exp(1i*dataPointsLine(4)*(pi/180)); %#ok<AGROW>
                        sch(idxNetParams, ports(row), ports(col_3)) = 10^(dataPointsLine(5)/20) * exp(1i*dataPointsLine(6)*(pi/180)); %#ok<AGROW>
                        sch(idxNetParams, ports(row), ports(col_4)) = 10^(dataPointsLine(7)/20) * exp(1i*dataPointsLine(8)*(pi/180)); %#ok<AGROW>
                    otherwise
                end
            end
            idxLine = idxLine+1;
            if idxLine>length(data_file)
                break;
            end
        end

        freqCounter = freqCounter+1; % capturing number of freq points for preallocation.
        idxNetParams = idxNetParams + 1;
    end
    if idxLine > length(data_file)
        break
    end
end

D=NaN(size(sch));
% calculate differential s parameter matrix from single ended
for i=1:size(sch,1)
    S(:,:) = sch(i,:,:);
    T = [1 1 0 0 ; 1 -1 0 0 ; 0 0 1 1 ; 0 0 1 -1];
    W = T * (S / T);
    D(i,:,:) = W(:,:);
end

% D matrix should be
% Scc11 Scd11 Scc12 Scd21
% Sdc11 Sdd11 Sdc12 Sdd12
% Scc21 Scd21 Scc22 Scd22
% Sdc21 Sdd21 Sdc22 Sdd22

% proper values
SDD(:,1,1) = D(:,2,2);
SDD(:,2,2) = D(:,4,4);
SDD(:,1,2) = D(:,2,4);
SDD(:,2,1) = D(:,4,2);

SDC(:,1,1) = D(:,2,1);
SDC(:,2,2) = D(:,4,3);
SDC(:,1,2) = D(:,2,3);
SDC(:,2,1) = D(:,4,1);

% backwards compatibility output variables
data.m = sch;
schFreqAxis=schFreqAxis(1:freqCounter); % truncating preallocated array to number of freq points.
data.freq = schFreqAxis;
colors = 'rgbk';

if (plot_ini_s_params == 1)
    figure('name', 'Single-ended s-parameters');
    for mj=1:4
        subplot(2,2,mj);
        for mi=1:4
            plot(data.freq, 20*log10(abs(data.m(:,mj,mi)+1.0e-15)), ...
                colors(mi), 'linewidth', 2, 'disp', sprintf('S%d%d', mj, mi));
            hold on
        end
        xlabel('Frequency (Hz)');
        ylabel('Magnitude (dB)');
        legend show
        grid on
        title(sprintf('Output port %d', mj));
    end
end

if (plot_dif_s_params == 1)
    figure('name', 'Mixed-mode s-parameters');
    subplot(2,1,1);
    for mj=1:2
        for mi=1:2
            plot(data.freq, 20*log10(abs(SDD(:,mj,mi))), ...
                colors((mj-1)*2+mi), 'linewidth',2, 'disp', sprintf('SDD%d%d', mj, mi));
            hold on
        end
    end
    xlabel('Frequency (Hz)');
    ylabel('Magnitude (dB)');
    legend show
    grid on
    title(infile);

    subplot(2,1,2);
    for mj=1:2
        for mi=1:2
            plot(data.freq, 20*log10(abs(SDC(:,mj,mi))+1.0e-15), ...
                colors((mj-1)*2+mi), 'linewidth',2, 'disp', sprintf('SDC%d%d', mj, mi));
            hold on
        end
    end
    xlabel('Frequency (Hz)');
    ylabel('Magnitude (dB)');
    legend show
    grid on
end

if OP.DISPLAY_WINDOW, close(hMsgBox); end
% end read_sp4_sparam
function [data, SDD, SDC] = read_p2_s2params(infile, plot_ini_s_params, plot_dif_s_params, ports,OP)
%% FUNCTION :: read_sp4_sparams
%
% Description
%   Read the fid of single-ended 4-port complex S-parameters
%   in Touchstone format 'file' and convert to the internal
%   format using the port transform 'ports'
%
%   Created by Mike Y. He
%   April 22, 2005
%
%   Reused some code from
%   Anthony Sanders, Alex Deas, Bob Davidov (24 January 2005)
%   for touchstone 4-port S-matrix import.
%
%   Modified (2012-July-27) by Ken Young to match current indexing scheme and
%   optimized for quicker parameter matching and parsing. also, separated out
%   the plotting algorithms into their own sub-function routines
%
% Input Variables (required)
%   file                -- The s4p file to be read and converted
%   plot_ini_sparams    -- Plot the initial s-parameter information. For debugging purposes
%   plot_diff_sparams   -- Plot the differential s-parameter information. For debugging purposes
%   swap                -- Re-order the port layout
%
% Output/Return Variables
%   sch         -- the data matrix contain the network parameter data points
%   schFreq     -- the frequency vector of the network parameter data points
%   sdc         -- the differential in/common-mode out s-parameter data matrix
%   sdd         -- the differential in/differential out s-parameter data matrix
%
% Local Variables (scoped locally to current function only)
%   fid                 -- the opened file being processed
%   fileDataLine        -- the current line in the fid
%   optionLine          -- the options line in the s4p file
%   ports               -- the order of the ports being processed
%   idxNetParams        -- the position (line) of the network parameter data points in the s4p file
%   netParamDataPoints  -- the network parameter data points of the current line of the fid
%   schFreq             -- the frequency points in the s4p file stored as a separate vector
%   sch                 -- the network parameter data points from the s4p file converted in matrix form
%   S                   --
%   T                   --
%   W                   --
%   D                   --
%   sdd                 -- the differential in/differential out form of network parameter data points
%   sdc                 -- the differential in/common-mode out form of network parameter data points
%   plot_ini_sparams    -- Plot the initial s-parameter information. For debugging purposes
%   plot_diff_sparams   -- Plot the differential s-parameter information. For debugging purposes
%

% open the fid
fid = fopen(infile, 'r');

data_file = textscan(fid, '%s', 'delimiter', sprintf('\n'));
data_file = data_file{1,1}; % removes a dimension layer of cell array depth

fclose(fid);

% backwards compatibility settings. can be removed in updated code.
if ~exist('OP', 'var'); OP.DISPLAY_WINDOW = true; end
if isempty(ports); ports = [1 2]; end % default order normally used.
ports = [1 2];


if OP.DISPLAY_WINDOW
    set(0,'defaulttextinterpreter','none') % prevents subscripting character in displayed messages
    hMsgBox = msgbox(infile, 'Reading S-Parameter File'); % display a progress bar for reading the s-parameter file(s)
end

% parse the file array until the s-parameters options line is found
idxLine = 1;
while true
    if regexp(data_file{idxLine}, '^#', 'once')
        % tests the option fileDataLine cell array to see if it contains...
        if ~isempty(regexpi(data_file{idxLine}, '\sGHz')); freq_scale = 1e9; end
        if ~isempty(regexpi(data_file{idxLine}, '\sMHz')); freq_scale = 1e6; end
        if ~isempty(regexpi(data_file{idxLine}, '\sKHz')); freq_scale = 1e3; end
        if ~isempty(regexpi(data_file{idxLine}, '\sHz'));  freq_scale = 1e0; end
        if ~isempty(regexpi(data_file{idxLine}, 'RI')) && ~isempty(regexpi(data_file{idxLine}, 'R'))
            format = 'complex';
        end
        if ~isempty(regexpi(data_file{idxLine}, 'MA')) && ~isempty(regexpi(data_file{idxLine}, 'R'))
            format = 'linpolar';
        end
        if ~isempty(regexpi(data_file{idxLine}, 'DB')) && ~isempty(regexpi(data_file{idxLine}, 'R'))
            format = 'dBpolar';
        end

        % if the format or freq_scale are not set then exit with error
        if isempty(freq_scale)
            error( 'read_s4p:MissingOptionsParameter', ['\n\t THE FREQUENCY SCALE PARAMETER INFORMATION IS MISSING' ...
                'FROM THE OPTIONS LINE IN THE S2P FILE']);
        elseif isempty(format)
            error( 'read_s4p:MissingOptionsParameter', ['\n\t THE FORMAT PARAMETER INFORMATION IS MISSING' ...
                'FROM THE OPTIONS LINE IN THE S2P FILE']);
        end

        idxLine = idxLine+1; %step one more line before exiting loop
        break
    end
    idxLine = idxLine+1;
end

col_1=1; col_2=2; % This is to maintain current indexing compatibility as well as readability
schFreqAxis = zeros(1,int8((length(data_file)-idxLine)/4)); % preallocated array length for memory and time purposes
idxNetParams = 1; % each 4-port network parameter data point set will receive an index

freqCounter = 0; % variable that will be used to capture number of freq points.

while true
    if isempty(data_file{idxLine}) || ~isempty(regexp(data_file{idxLine}, '^!', 'once'))
        idxLine = idxLine+1;
    else
        for row=1:2
%             dataPointsLine = eval (['[' data_file{idxLine} ']']);
            dataPointsLine= sscanf(data_file{idxLine}, '%f'); % Close V2015 memory leak - Adam Gregory
            if mod(length(dataPointsLine),4)==1
                schFreqAxis(idxNetParams) = dataPointsLine(1) * freq_scale;
                dataPointsLine(1) = []; % removes the frequency value from the current data fileDataLine cell array
                % so that only the network parameter data is in the array
            end
            if mod(length(dataPointsLine),4)~=0  % sanity check to make sure the network parameter data point line structure is correct
                msg = sprintf(['\n'...
                    '\n\t ******************************************************' ...
                    '\n\t ** There is an error in the s4p network parameter   **' ...
                    '\n\t ** data point information structure at line %G.'        ...
                    '\n\t **              ENDING PROGRAM                      **' ...
                    '\n\t ******************************************************' ...
                    '\n\n\n\t Press any key to exit (an error will be displayed). \n\n']);
                disp(msg)
                pause on;
                pause;
                pause off;
                error( 'COM:reads_4p:SparamFileDataError', 'INCORRECT S4P NETWORK PARAMETER DATA POINT STRUCTURE');
            else
                switch format
                    case 'complex'
                        sch(idxNetParams, ports(row), ports(col_1)) = dataPointsLine(1) + 1i*dataPointsLine(2); %#ok<AGROW>
                        sch(idxNetParams, ports(row), ports(col_2)) = dataPointsLine(3) + 1i*dataPointsLine(4); %#ok<AGROW>
                    case 'linpolar'
                        sch(idxNetParams, ports(row), ports(col_1)) = dataPointsLine(1) * exp(1i*dataPointsLine(2)*(pi/180)); %#ok<AGROW>
                        sch(idxNetParams, ports(row), ports(col_2)) = dataPointsLine(3) * exp(1i*dataPointsLine(4)*(pi/180)); %#ok<AGROW>
                    case 'dBpolar'
                        sch(idxNetParams, ports(row), ports(col_1)) = 10^(dataPointsLine(1)/20) * exp(1i*dataPointsLine(2)*(pi/180)); %#ok<AGROW>
                        sch(idxNetParams, ports(row), ports(col_2)) = 10^(dataPointsLine(3)/20) * exp(1i*dataPointsLine(4)*(pi/180)); %#ok<AGROW>
                    otherwise
                end
            end
            idxLine = idxLine+1;
            if idxLine>length(data_file)
                break;
            end
        end

        freqCounter = freqCounter+1; % capturing number of freq points for preallocation.
        idxNetParams = idxNetParams + 1;
    end
    if idxLine > length(data_file)
        break
    end
end

D=NaN(size(sch));
% calculate differential s parameter matrix from single ended
for i=1:size(sch,1)
    S(:,:) = sch(i,:,:);
    T = [1 1  ; 1 -1  ];
    W = T * (S / T);
    D(i,:,:) = W(:,:);
end

% D matrix should be
% Scc11 Scd11 Scc12 Scd21
% Sdc11 Sdd11 Sdc12 Sdd12
% Scc21 Scd21 Scc22 Scd22
% Sdc21 Sdd21 Sdc22 Sdd22

% proper values
SDD(:,1,1) = D(:,1,2);
SSC(:,1,1)=  D(:,2,1);



% backwards compatibility output variables
data.m = sch;
schFreqAxis=schFreqAxis(1:freqCounter); % truncating preallocated array to number of freq points.
data.freq = schFreqAxis;
colors = 'rgbk';

if (plot_ini_s_params == 1)
    figure('name', 'Single-ended s-parameters');
    for mj=1:4
        subplot(2,2,mj);
        for mi=1:4
            plot(data.freq, 20*log10(abs(data.m(:,mj,mi)+1.0e-15)), ...
                colors(mi), 'linewidth', 2, 'disp', sprintf('S%d%d', mj, mi));
            hold on
        end
        xlabel('Frequency (Hz)');
        ylabel('Magnitude (dB)');
        legend show
        grid on
        title(sprintf('Output port %d', mj));
    end
end

if (plot_dif_s_params == 1)
    figure('name', 'Mixed-mode s-parameters');
    subplot(2,1,1);
    for mj=1:2
        for mi=1:2
            plot(data.freq, 20*log10(abs(SDD(:,mj,mi))), ...
                colors((mj-1)*2+mi), 'linewidth',2, 'disp', sprintf('SDD%d%d', mj, mi));
            hold on
        end
    end
    xlabel('Frequency (Hz)');
    ylabel('Magnitude (dB)');
    legend show
    grid on
    title(infile);

    subplot(2,1,2);
    for mj=1:2
        for mi=1:2
            plot(data.freq, 20*log10(abs(SDC(:,mj,mi))+1.0e-15), ...
                colors((mj-1)*2+mi), 'linewidth',2, 'disp', sprintf('SDC%d%d', mj, mi));
            hold on
        end
    end
    xlabel('Frequency (Hz)');
    ylabel('Magnitude (dB)');
    legend show
    grid on
end

if OP.DISPLAY_WINDOW, close(hMsgBox); end
% end read_sp4_sparam

function result = readdataSnPx(filename, nport)
%function [freq, cs] = readdataSnPx(filename, nport)
% [freq, cs] = readdataSnP(filename, nport, format, nheader)
%
% Read Touchstone file with frequencies in units of Hertz
%
% Input:
% ======
% filename: Name of the Touchstone/SnP file
% nport: Number of ports
% format: 'RI' for real/imag, 'MA' for mag/angle (check option line in the
%         Touchstone file)
% nheader: Number of header lines (comment lines plus option line in the
%          Touchstone file)
%
% Output:
% =======
% freq: Vector of frequencies [Hz]
% cs: 3-D array of complex-valued S parameters where cs(i,j,k) is S(i,j)
%     at frequency freq(k)
%
% Note: If frequency unit is not Hertz (but GHz, MHz etc.) simply scale
% frequencies appropriately after reading the data.
%
% Ref.: Touchstone(R) File Format Specification, Rev.1.1,
% EIA/IBIS Open Forum, 2002.
%
% Written by Henning Braunisch, September 2004.
% Updated by Steven Krooswyk, April 2006.


fid = fopen(filename, 'r');


% Skip header lines
str = ' ';
n = 0;
while ~strcmp(str(1),'#')
    str = fgetl(fid);
    if isempty(str)
        str=' ' ;
        if n > 1000
            display('error: could not find config line (#)')
            break
        end
    end
    n = n + 1;
end

% parse configuration line
A=sscanf(str,'%1s %2s %1s %2s %1s %2s',[1,inf]);
p = find(A=='S');           %position of 'S'
units = lower(A(2:p-1));    %units before 'S'
format = A(p+1:p+2);        %format after 'S'

% skip any more header lines
%while ~str

nk = 0; % frequency counter
while 1

    [temp, count] = fscanf(fid, '%f', 1);
    if count == 0
        temp2 = fscanf(fid, '%s', 1);
        if ~isempty(temp2), fgetl(fid); continue, end;
        break
    end
    nk = nk+1; freq(1,nk) = temp;  %#ok<AGROW>
    for ni = 1:nport
        for nj = 1:nport
            switch lower(format)
                case 'ma'
                    mag = fscanf(fid, '%f', 1); ang = fscanf(fid, '%f', 1);
                    cs(ni,nj,nk) = mag * exp(1i*ang*pi/180); %#ok<AGROW>
                case 'ri'
                    re = fscanf(fid, '%f', 1); im = fscanf(fid, '%f', 1);
                    cs(ni,nj,nk) = complex(re, im); %#ok<AGROW>
                case 'db'
                    db = fscanf(fid, '%f', 1); ang = fscanf(fid, '%f', 1);
                    M = 10^(db/20);
                    %re = M*cos(ang);
                    %im = M*sin(ang);
                    re =  M*cos(ang * pi / 180);
                    im =  M*sin(ang * pi / 180);
                    cs(ni,nj,nk) = complex(re, im); %#ok<AGROW>
                otherwise
                    error('readdataSnP: Unknown data format');
            end
        end
    end
end

fclose(fid);

% If 2-port then swap S_12 and S_21 per Touchstone spec
if nport == 2
    temp = cs(2,1,:);
    cs(2,1,:) = cs(1,2,:);
    cs(1,2,:) = temp;
end

% Update freq units to Hz
switch lower(units)
    case 'hz'

    case 'khz'
        freq=freq.*1e3;
    case 'mhz'
        freq=freq.*1e6;
    case 'ghz'
        freq=freq.*1e9;
end

% passivity check
result.freq = freq;
result.cs    = cs;

function result = reduce(var1)
% --- Reduce 1x1xn array to 1xn (aka squeeze)
out = zeros(1,length(var1));
out(1,:) = var1(1,1,:);
result=out;

function pdf=get_pdf(chdata, delta_y, t_s, param, OP)
SBR=chdata.eq_pulse_response(:)'; % row vector
type=chdata.type;
samp_UI=param.samples_per_ui;
residual_response = SBR;

if isequal(type, 'THRU')
    % for thru pulse response:
    % remove the cursor and the DFE postcursors (up to their limit), since
    % we only care about the residuals.
    ideal_cancelled_cursors = SBR(t_s+param.samples_per_ui*(0:param.ndfe));
    effective_cancelled_cursors = sign(ideal_cancelled_cursors) .* ...
        min(abs(ideal_cancelled_cursors), residual_response(t_s)*[1, param.bmax]);
    effective_cancellation_samples = kron(effective_cancelled_cursors, ones(1, param.samples_per_ui));

    % Apply a constant DFE coefficient 1/2 UI before and after each postcursor. Not
    % really needed for COM, but helps debugging. May be factored out in future revisions.
    start_cancel = t_s-param.samples_per_ui/2;
    end_cancel = t_s+(1/2+param.ndfe)*param.samples_per_ui - 1;
    
    residual_response(start_cancel:end_cancel) = ...
        residual_response(start_cancel:end_cancel) - effective_cancellation_samples;
    %else
    % for crosstalk pulse responses, nothing is cancelled, and all phases
    % are equally important.
end

nui=round(length(residual_response)/param.samples_per_ui);

vs=zeros(nui-2, param.samples_per_ui);
for i=1:param.samples_per_ui
    vs(:,i)=residual_response(param.samples_per_ui*(1:nui-2)+i);
end

if OP.DISPLAY_WINDOW,
    hwaitbar=waitbar(0);
end

% determine which pdf to use
if isequal(type, 'THRU')
    % one phase is interesting for thru
    phases = mod(t_s,param.samples_per_ui);
    if phases==0, phases = param.samples_per_ui; end
else
    phases=1:samp_UI;
end

mxV = zeros(size(phases));
for k=phases
    pdf_samples(k)=get_pdf_from_sampled_signal(vs(:,k), param.levels, delta_y); %#ok<AGROW>
    mxV(k)=sqrt(sum( pdf_samples(k).x.^2.*pdf_samples(k).y)); % standard deviation of PDF
    progress = k/length(phases);
    if OP.DISPLAY_WINDOW, waitbar(progress, hwaitbar, ['processing COM pdf ' chdata.base ] ); figure(hwaitbar); drawnow; end
end
[UNUSED_OUTPUT pxi]=max(mxV); %#ok<ASGLU>
pdf=pdf_samples(pxi);

if OP.DISPLAY_WINDOW,
    close(hwaitbar);
end

function [ILN, efit]= get_ILN(sdd21,faxis_f2)
% used for FD IL fitting
% sdd21 us a complex insertion loss
db = @(x) 20*log10(abs(x));
sdd21=squeeze(sdd21);
if  iscolumn(sdd21)
    sdd21=sdd21.';
end
fmbg=[ones(length(faxis_f2),1).*transpose(abs(sdd21))  transpose(sqrt(faxis_f2)).*transpose(abs(sdd21))  transpose(faxis_f2).*transpose(abs(sdd21)) transpose(faxis_f2.^2).*transpose(abs(sdd21)) ];
warning('off','MATLAB:nearlySingularMatrix');
LGw=transpose(abs(sdd21).*db(sdd21));
alpha = ((fmbg'*fmbg)^-1)*fmbg'*LGw;
efit=(alpha(1)+alpha(2).*sqrt(faxis_f2)+alpha(3).*faxis_f2 +faxis_f2.^2.*alpha(4)   );
ILN = db(sdd21)-efit;

function result=optimize_fom(OP, param, chdata, sigma_bn)
%% input
% chdata(1).uneq_imp_response is the impulse response input expected to be normalized to At, peak drive voltage
% baud_rate - baud rate in seconds
% param.samples_per_ui = samples per UI of IR
% param.max_ctle - maximum ac to dc gain in dB
% param.tx_ffe(1) - maximum pre cursor (positive value)
% param.tx_ffe(2) - maximum post cursor (positive value)
% param.tx_ffe_step - sweep step size for tx pre and post taps
% param.ndfe - number of reference dfe taps
% output
% result.eq.txle - [ precusor curosr postcursor]: pre and post are negative
% result.eq.ctle - index of CTLE parameters in table
% result.IR - impulse response
% result.avail_signal - maximum signal after equalization
% result.avail_sig_index - index in result.IR of max signal
% result.best_FOM - best raw ISI

min_number_of_UI_in_response=40;
baud_rate=1/param.ui;

H_r = 1./polyval([1 2.613126 3.414214 2.613126 1], 1i*chdata(1).faxis/(param.f_r*param.fb));
cm4_values = param.tx_ffe_cm4_values;
cm3_values = param.tx_ffe_cm3_values;
cm2_values = param.tx_ffe_cm2_values;
cm1_values = param.tx_ffe_cm1_values;
cp1_values = param.tx_ffe_cp1_values;

gdc_values = param.ctle_gdc_values;
switch param.CTLE_type
    case 'CL93'
    case 'CL120d'
        g_DC_HP_values =param.g_DC_HP_values;
    case 'CL120e'
        f_HP_Z=param.f_HP_Z;
        f_HP_P=param.f_HP_P;
        
end
best_ctle = [];
best_FOM = -inf;
best_txffe = [];
delta_sbr = [];
pxi=0;
if OP.DISPLAY_WINDOW
    hwaitbar=waitbar(0);
else
    fprintf('FOM search');
end
FOM=0;
for ctle_index=1:length(gdc_values)
    g_dc = gdc_values(ctle_index);
    kacdc = 10^(g_dc/20);
    CTLE_fp1 = param.CTLE_fp1(ctle_index);
    CTLE_fp2 = param.CTLE_fp2(ctle_index);
    CTLE_fz = param.CTLE_fz(ctle_index);
    switch param.CTLE_type
        case 'CL93'
            lf_indx=1;
        case 'CL120d'
            lf_indx=length(g_DC_HP_values);
        case 'CL120e'
            lf_indx=1;
            HP_Z = param.f_HP_Z(ctle_index);
            HP_P = param.f_HP_P(ctle_index);
    end
    for  g_LP_index=1:lf_indx
        ctle_gain = (kacdc + 1i*chdata(1).faxis/CTLE_fz) ./ ...
            ((1+1i*chdata(1).faxis/CTLE_fp1).*(1+1i*chdata(1).faxis/CTLE_fp2));
        switch param.CTLE_type
            case 'CL93'
                H_low=1;
                kacde_DC_low=1;
            case 'CL120d'
                g_DC_low = g_DC_HP_values(g_LP_index);
                f_HP=param.f_HP(g_LP_index);
                kacde_DC_low = 10^(g_DC_low/20);
                H_low=(kacde_DC_low +  1i*chdata(1).faxis/f_HP)./(1 + 1i*chdata(1).faxis/f_HP);
            case 'CL120e' % z1 has been adusted on read in
                H_low=(1 +  1i*chdata(1).faxis/HP_Z)./(1 + 1i*chdata(1).faxis/HP_P);
        end
        H_ctf=H_low.*ctle_gain;
        if OP.INCLUDE_CTLE==1
            for k=1:param.num_s4p_files
                chdata(k).ctle_imp_response = TD_CTLE(chdata(k).uneq_imp_response, baud_rate ...
                    , CTLE_fz, CTLE_fp1, CTLE_fp2, g_dc, param.samples_per_ui);
                switch param.CTLE_type
                    case 'CL93'
                    case 'CL120d'
                        chdata(k).ctle_imp_response = TD_CTLE(chdata(k).ctle_imp_response, baud_rate, f_HP, f_HP,100e100 , g_DC_low , param.samples_per_ui);
                    case 'CL120e' % z1 has been adusted on read in
                         chdata(k).ctle_imp_response = TD_CTLE(chdata(k).ctle_imp_response, baud_rate, HP_Z,HP_P,100e100 , 0 , param.samples_per_ui);
                end
            end
 
        else
            for k=1:param.num_s4p_files
                chdata(k).ctle_imp_response = chdata(k).uneq_imp_response;
            end
        end
        
        %% Equation 93A-22 %%
%         figure(1000)
%         semilogx(chdata(1).faxis/1e9,db(H_ctf))
%         hold on
        if OP.RX_CALIBRATION
            ctle_gain2 = (kacdc + 1i*chdata(2).faxis/CTLE_fz) ./ ...
                ((1+1i*chdata(2).faxis/CTLE_fp1).*(1+1i*chdata(2).faxis/CTLE_fp2));
            
            
        switch param.CTLE_type
            case 'CL93'
                H_low2=1;
            case 'CL120d'
                g_DC_low = g_DC_HP_values(g_LP_index);
                f_HP=param.f_HP(g_LP_index);
                kacde_DC_low = 10^(g_DC_low/20);
                H_low2=(kacde_DC_low +  1i*chdata(1).faxis/f_HP)./(1 + 1i*chdata(1).faxis/f_HP);
            case 'CL120e' % z1 has been adusted on read in
                H_low2=(1 +  1i*chdata(1).faxis/HP_Z)./(1 + 1i*chdata(1).faxis/HP_P);
        end
            H_ctf2=H_low2.*ctle_gain2;
        end
        %% Equation 93A-35 - independent of FFE setting %%
        sigma_N = sqrt(param.eta_0*sum(abs(H_r(2:end) .* H_ctf(2:end)).^2 .* diff(chdata(1).faxis)/1e9));
        if OP.RX_CALIBRATION
            sigma_ne = get_sigma_noise( H_ctf2, param, chdata, sigma_bn); %% Equation 93A-48 %%
            sigma_NEXT=0;
        else
            %% Equations 93A-33 and 93A-34  for NEXT - independent of TXFFE setting %%
            sigma_NEXT =  get_xtlk_noise( [0 1 0], 'NEXT', param, chdata );
            sigma_ne=0;
        end
        for k_cm4=1:length(cm4_values)
            cm4=cm4_values(k_cm4);
            for k_cm3=1:length(cm3_values)
                cm3=cm3_values(k_cm3);
                for k_cm2=1:length(cm2_values)
                    cm2=cm2_values(k_cm2);
                    for k_cm1=1:length(cm1_values)
                        cm1=cm1_values(k_cm1);
                        for k_cp1=1:length(cp1_values)
                            cp1=cp1_values(k_cp1);
                            pxi=pxi+1;
                            progress = pxi/(length(cm2_values)* length(cm1_values)*length(cp1_values)*length(gdc_values)*length(cm3_values)*length(cm4_values)*lf_indx );
                            txffe = [cm4, cm3, cm2, cm1, 1-abs(cm3)-abs(cm4)-abs(cm2)-abs(cm1)-abs(cp1), cp1];
                            % Skip combinations with small values of c(0), not guaranteed to be supported by all transmitters.
                            if txffe(5)<param.tx_ffe_c0_min
                                continue;
                            end
                            upsampled_txffe = zeros(1, param.samples_per_ui*5+1); % start with zeros everywhere
                            upsampled_txffe(1+(0:5)*param.samples_per_ui) = txffe; % "plant" the coefficients in the desired locations
                            effective_channel = filter(upsampled_txffe, 1, chdata(1).ctle_imp_response);
                            sbr=filter(ones(param.samples_per_ui, 1), 1, effective_channel);
                            [UNUSED_OUTPUT, sbr_peak_i]=max(abs(sbr)); %#ok<ASGLU>
                            triple_transit_time = round(sbr_peak_i*2/param.samples_per_ui)+20;
                            if min_number_of_UI_in_response < triple_transit_time
                                min_number_of_UI_in_response = triple_transit_time;
                            end
                            
                            %% initial guess at cursor location (t_s)  - based on approximate zero crossing
                            zxi = find(diff(sign(sbr-.01*max(sbr)))>=1);
                            zxi = zxi(zxi<sbr_peak_i);
                            zxi = zxi(sbr_peak_i - zxi < 4*param.samples_per_ui);
                            if isempty(zxi)
                                continue;
                            elseif length(zxi)>1
                                zxi=zxi(end);
                            end
                            
                            if param.ndfe==0
                                max_dfe1=0;
                            else
                                max_dfe1=param.bmax(1);
                            end
                            %% adjust cursor_i to Solve equation 93A-25 %%
                            % Muller-Mueller criterion with DFE
                            mm_range = zxi+(0:2*param.samples_per_ui);
                            mm_metric = ...
                                abs(sbr(mm_range-param.samples_per_ui) - max(sbr(mm_range+param.samples_per_ui)-max_dfe1*sbr(mm_range), 0));
                            [UNUSED_OUTPUT, mm_cursor_offset] = min(mm_metric); %#ok<ASGLU>
                            cursor_i = zxi+mm_cursor_offset-1;
                            
                            cursor = sbr(cursor_i);
                            %% 93A.1.6 step c defines A_s %%
                            A_s = param.R_LM*cursor/(param.levels-1);
                            if isempty(delta_sbr)
                                delta_sbr = sbr;
                            end
                            sbr=sbr(:);
                            
                            %% Equation 93A-27 "otherwise" case %%
                            far_cursors = sbr(cursor_i+param.samples_per_ui*(param.ndfe+1):param.samples_per_ui:end);
                            t=((cursor_i+param.samples_per_ui*(param.ndfe+1):param.samples_per_ui:length(sbr))-(cursor_i+param.samples_per_ui*(param.ndfe+1)))*...
                                param.ui/param.samples_per_ui;
                            precursors = sbr(cursor_i-param.samples_per_ui:-param.samples_per_ui:1);
                            precursors = precursors(end:-1:1);
                            
                            % Error message if the sbr is not long enough for the specified range of Nb
                            if length(sbr) < cursor_i+param.samples_per_ui*(param.ndfe+1)
                                close(hwaitbar);
                                error('Pulse Response contains %d samples after the cursor. Specified Nb requires %d samples after the cursor.' ...
                                    , length(sbr)-cursor_i, param.samples_per_ui*(param.ndfe+1));
                            end
                            
                            %% Equation 93A-27, when 1<=n<=N_b
                            dfecursors=sbr(cursor_i+param.samples_per_ui*(1):param.samples_per_ui:cursor_i+param.samples_per_ui*(param.ndfe));
                            excess_dfe_cursors = dfecursors - ...
                                sign(dfecursors) .*min(sbr(cursor_i)*param.bmax(:), abs(dfecursors) );
                            
                            %% Eq. 93A-28 %%
                            sampling_offset = mod(cursor_i, param.samples_per_ui);
                            %ensure we can take early sample
                            if sampling_offset<=1
                                sampling_offset=sampling_offset+param.samples_per_ui;
                            end
                            if (OP.LIMIT_JITTER_CONTRIB_TO_DFE_SPAN)
                                cursors_early_sample = sbr(cursor_i-1+param.samples_per_ui*(-1:param.ndfe));
                                cursors_late_sample = sbr(cursor_i+1+param.samples_per_ui*(-1:param.ndfe));
                            else
                                cursors_early_sample = sbr(sampling_offset-1:param.samples_per_ui:end);
                                cursors_late_sample = sbr(sampling_offset+1:param.samples_per_ui:end);
                            end
                            % ensure lengths are equal
                            cursors_early_sample = cursors_early_sample(1:length(cursors_late_sample));
                            h_J = (cursors_late_sample-cursors_early_sample)/2*param.samples_per_ui;
                            
                            %% Equation 93A-30 %%
                            % since A_s = param.R_LM*cursor/(param.levels-1), cursor=(param.levels-1)*A_s/param.R_LM
                            sigma_TX = (param.levels-1)*A_s/param.R_LM*10^(-param.SNR_TX/20);
                            
                            %% Equation 93A-31 %%
                            sigma_ISI = param.sigma_X*norm([precursors; excess_dfe_cursors; far_cursors]);
                            
                            %% Equation 93A-32 %%
                            sigma_J = norm([param.A_DD param.sigma_RJ])*param.sigma_X*norm(h_J);
                            
                            %% Equations 93A-33 and 93A-34 for FEXT (depends on TXFFE setting) %%
                            if OP.RX_CALIBRATION
                                sigma_XT=0;
                            else
                                sigma_FEXT =  get_xtlk_noise( upsampled_txffe, 'FEXT', param, chdata );
                                sigma_XT = norm([sigma_NEXT sigma_FEXT]);
                            end
                            %% Equation 93A-36 denominator (actually its sqrt)
                            total_noise_rms = norm([sigma_ISI sigma_J sigma_XT sigma_N sigma_TX sigma_ne]);
                            
                            %% Equation 93A-36 (note log argument is voltage rather than power ratio)
                            FOM = 20*log10(A_s/total_noise_rms);
                            
                            if (FOM > best_FOM)
                                best_txffe = txffe;
                                best_sbr = sbr;
                                best_ctle = ctle_index;
                                best_G_high_pass =g_LP_index;
                                best_FOM = FOM;
                                best_cursor_i = cursor_i;
                                best_IR=effective_channel;
                                best_sigma_N = sigma_N;
                                best_h_J = h_J;
                                best_A_s=A_s;
                            end
                        end
                        if OP.DISPLAY_WINDOW
                            waitbar(progress, hwaitbar, 'Linear equalization tuning'); figure(hwaitbar); drawnow;
                        else
                            if ~mod(progress*100,20), fprintf('%i%% ', progress*100 );end
                        end
                    end
                end
            end
        end
    end
end



if ~exist('best_cursor_i', 'var')% take last setting
    best_txffe = txffe;
    best_sbr = sbr;
    best_ctle = ctle_index;
    best_G_high_pass =g_LP_index;
    best_FOM = FOM;
    best_cursor_i = cursor_i;
    best_IR=effective_channel;
    best_sigma_N = sigma_N;
    best_h_J = h_J;
end

best_cursor = best_sbr(best_cursor_i);
% report during debug
PRin=filter(ones(param.samples_per_ui, 1),1, chdata(1).uneq_imp_response);
f=1e8:1e8:100e9;
H_r = 1./polyval([1 2.613126 3.414214 2.613126 1], 1i*f./(param.f_r*param.fb));
ctle_gain1 = (10^(gdc_values(best_ctle)/20) + 1i*f/param.CTLE_fz(best_ctle)) ./ ...
    ((1+1i*f/param.CTLE_fp1(best_ctle)).*(1+1i*f/param.CTLE_fp2(best_ctle)));

switch param.CTLE_type
    case 'CL93'
        H_low=1;
    case 'CL120d'
        H_low=(10^(param.g_DC_HP_values(best_G_high_pass)/20) +  1i*f/param.f_HP(best_G_high_pass))./(1 + 1i*f/param.f_HP(best_G_high_pass));
    case 'CL120e'
        H_low=(1 +  1i*f/param.f_HP_Z(best_ctle))./ (1 + 1i*f/param.f_HP_P(best_ctle));
end
ctle_gain=H_low.*ctle_gain1.*H_r;



lsbr=length(sbr);
t=0:param.ui/param.samples_per_ui:(lsbr-1)*param.ui/param.samples_per_ui;

sampled_best_sbr_precursors_t = (best_cursor_i/param.samples_per_ui:-1:1/param.samples_per_ui)*param.ui;
sampled_best_sbr_precursors_t = sampled_best_sbr_precursors_t(end:-1:2); % exclude cursor
sampled_best_sbr_precursors = best_sbr(round(sampled_best_sbr_precursors_t/param.ui*param.samples_per_ui));
sampled_best_sbr_postcursors_t = (best_cursor_i:param.samples_per_ui:lsbr)/param.samples_per_ui*param.ui;
sampled_best_sbr_postcursors_t = sampled_best_sbr_postcursors_t(2:end); % exclude cursor
sampled_best_sbr_postcursors = best_sbr(round(sampled_best_sbr_postcursors_t/param.ui*param.samples_per_ui));
sampled_best_sbr_dfecursors_t = (best_cursor_i/param.samples_per_ui+(1:param.ndfe))*param.ui;
% apply max tap value constraint
dfe_cursors = sampled_best_sbr_postcursors(1:param.ndfe);
DFE_taps_mV = sign(dfe_cursors).*min(best_cursor*param.bmax(:), abs(dfe_cursors) );
sampled_best_sbr_postcursors(1:param.ndfe) = dfe_cursors-DFE_taps_mV;

if OP.DEBUG ~=0
    display(['FOM:                ' ,num2str(best_FOM, 2),' dB']);
    display(['TXFFE coefficients: ' ,mat2str(best_txffe) ] );
    display(['CTLE DC gain:       ' ,num2str(gdc_values(best_ctle)), ' dB']);
    display(['CTF peaking gain:  ' ,num2str(20*log10(max(abs(ctle_gain))), 2), ' dB']);
    display(['Available signal:   ' ,num2str(best_cursor)]);
    
    if OP.DISPLAY_WINDOW && ~OP.RX_CALIBRATION
        figure(110);
        set(gcf, 'Name', 'CTLE selection');
        movegui(gcf, 'southeast');
        
        semilogx(f,20*log10(abs(ctle_gain)), 'disp', sprintf('Case %d', param.package_testcase_i));
        hold on
        fbaud_tick=find(f >= baud_rate, 1);
        fnq_tick=find(f >= baud_rate/2, 1);
        stem(f(fnq_tick),20*log10(abs(ctle_gain(fnq_tick))),'g', 'handlevisibility', 'off');
        stem(f(fbaud_tick),20*log10(abs(ctle_gain(fbaud_tick))),'g', 'handlevisibility', 'off');
        recolor_plots(gca);
        title('CTF/w Rx Filter Response')
        ylabel('dB')
        xlabel('Hz')

        % display pulse responses in one axis per test case.
        figure_name = sprintf('Equalization effect: %s', param.base);
        fig=findobj('Name', figure_name);
        if isempty(fig), fig=figure('Name', figure_name); end
        figure(fig);
        movegui(fig,'north')
        %figure(fig.Number);
        hax = subplot(length(OP.pkg_len_select), 1, param.package_testcase_i);

        plot(t,best_sbr,'disp', 'Equalized end-to-end PR');
        hold on
        PRplt(1:param.samples_per_ui+5)=PRin(1); % line up with the ffe introduced delay
        PRplt(param.samples_per_ui+6:length(t))=PRin(1:length(t)-param.samples_per_ui-5);
        plot(t,PRplt,'r','disp', 'Unequalized end-to-end PR');
        stem(t(best_cursor_i),best_sbr(best_cursor_i),'g','disp','Cursor (sample point)');
        title(sprintf('Case %d', param.package_testcase_i));
        ylabel('volts')
        xlabel('seconds')

        plot(sampled_best_sbr_precursors_t-param.ui/param.samples_per_ui, sampled_best_sbr_precursors', 'kx', 'disp','Pre cursors');
        plot(sampled_best_sbr_postcursors_t-param.ui/param.samples_per_ui, sampled_best_sbr_postcursors', 'ko', 'disp','Post cursors');
        stem(sampled_best_sbr_dfecursors_t-param.ui/param.samples_per_ui,dfe_cursors','m', 'LineWidth',2,'disp','DFE-canceled cursors');
        grid on
        legend show
        zoom xon
        set(hax, 'tag', 'EQE');
    end
end
if OP.DISPLAY_WINDOW
    close(hwaitbar);
else
    fprintf('\n');
end

% % eq_data
result.txffe = best_txffe;
result.ctle = best_ctle;
result.best_G_high_pass=best_G_high_pass;
result.DFE_taps = DFE_taps_mV/best_cursor; %relative
result.DFE_taps_i = best_cursor_i+(1:param.ndfe)*param.samples_per_ui;

result.IR = best_IR;
result.A_s = best_A_s;

result.t_s = best_cursor_i;
result.sigma_N = best_sigma_N;
result.h_J = best_h_J;
result.FOM = best_FOM;

function [impulse_response, p1_ctle, p2_ctle, z_ctle] = TD_CTLE(ir_in, fb, f_z, f_p1, f_p2, kacdc_dB, oversampling)
%% Equation 93A-22 implemented in z-space and applied to the impulse response.
p1_ctle = -2*pi*f_p1;
p2_ctle = -2*pi*f_p2;
z_ctle = -2*pi*f_z*10^(kacdc_dB/20);
k_ctle = -p2_ctle;
bilinear_fs = 2*fb*oversampling;
p2d = (1+p2_ctle/bilinear_fs)./(1-p2_ctle/bilinear_fs);
p1d = (1+p1_ctle/bilinear_fs)./(1-p1_ctle/bilinear_fs);
zd = (1+z_ctle/bilinear_fs)./(1-z_ctle/bilinear_fs);
% kd = (bilinear_fs-z_ctle)/((bilinear_fs-p1_ctle)*(bilinear_fs-p2_ctle));
% allow for different pole zeros RIM 9-29-2015
kd = (bilinear_fs-z_ctle)./((bilinear_fs-p1_ctle)*(bilinear_fs-p2_ctle))*f_p1/f_z;
B_filt =k_ctle*kd*poly([zd, -1]);
A_filt=poly([p1d, p2d]);
impulse_response=filter(B_filt,A_filt,ir_in);

function [ pdf ] = get_pdf_from_sampled_signal( input_vector, L, BinSize ,FAST_NOISE_CONV)
% Create PDF from interference vector using successive delta-set convolutions.
%   input_vector = list of values of samples
%   return
%   pdf.x
%   pdf.y
%   pdf.vec
%   pdf.bin
if ~exist('FAST_NOISE_CONV','var')
    FAST_NOISE_CONV=0;
end
if max(input_vector) > BinSize
    input_vector=input_vector(abs(input_vector)>BinSize);
end
% for i = 1:length(input_vector)
%    if abs(input_vector(i)) < BinSize , input_vector(i)=0; end
%end

input_vector(abs(input_vector)<BinSize) = 0;
b=sign(input_vector);
[input_vector,index]=sort(abs(input_vector),'descend');
input_vector=input_vector.*b(index);
if FAST_NOISE_CONV
    sig_res=norm(input_vector(find(abs(input_vector)<.001,1)+1:end));
    res_pdf= normal_dist(sig_res,5,BinSize);
    input_vector=input_vector(1:find(abs(input_vector)<.001,1));
end
%% Equation 93A-39 %%
values = 2*(0:L-1)/(L-1)-1;
prob = ones(1,L)/L;

%% Initialize pdf to delta at 0
pdf=d_cpdf(BinSize, 0, 1);

for k = 1:length(input_vector)
    pdfn=d_cpdf(BinSize, abs(input_vector(k))*values, prob);
    pdf=conv_fct(pdf, pdfn);
end
if FAST_NOISE_CONV
    pdf=conv_fct(pdf,res_pdf);
end

function pdf=d_cpdf( binsize, values, probs)
%  p=cpdf(type, ...)
%
% CPDF is a probability mass function for discrete distributions or an 
% approxmation of a PDF for continuous distributions.
% 
% cpdf is internally normalized so that the sum of probabilities is 1
% (regardless of bin size).

% Internal fields:
% Min: *bin number* of minimum value.
% BinSize: size of PDF bins. Bin center is the representative value.
% Vec: vector of probabilities per bin.

values=binsize*round(values/binsize);
t=(min(values):binsize:max(values));
pdf.Min=min(values)/binsize;
pdf.y=zeros(size(t));
for k=1:length(values)
    [UNUSED_OUTPUT, bin]=min(abs(t-values(k))); %#ok<ASGLU>
    pdf.y(bin) = pdf.y(bin)+probs(k);
end

pdf.BinSize=binsize;
pdf.y=pdf.y/sum(pdf.y);
if any(~isreal(pdf.y)) || any(pdf.y<0)
    error('PDF must be real and nonnegative');
end
support=find(pdf.y);
pdf.y=pdf.y(support(1):support(end));
pdf.Min=pdf.Min+(support(1)-1);
pdf.x=(pdf.Min:-pdf.Min)*binsize;

function p=conv_fct(p1, p2)
if p1.BinSize ~= p2.BinSize
    error('bin size must be equal')
end

p=p1;
p.BinSize=p1.BinSize;
%p.Min=p1.Min+p2.Min;
p.Min=round(p1.Min+p2.Min);	% modified by Yasuo Hidaka, 9/4/2016
p.y=conv(p1.y, p2.y);
%p.x =p.Min*p.BinSize:p.BinSize:-p.Min*p.BinSize;
p.x =(p.Min:-p.Min)*p.BinSize;	% modified by Yasuo Hidaka, 9/4/2016



function pdf = normal_dist(sigma,nsigma,binsize)
pdf.BinSize=binsize;
pdf.Min=-round(nsigma*sigma/binsize);
pdf.x=(pdf.Min:-pdf.Min)*binsize;
pdf.y=exp(-pdf.x.^2/(2*sigma^2+eps));
pdf.y=pdf.y/sum(pdf.y);

function [ param OP ]= read_ParamConfigFile(paramFile)
%warning('off','MATLAB:xlsread:Mode'); % suppress warning messages for reading the settings file from XLS
try
    [status,sheets] = xlsfinfo(paramFile); %#ok<ASGLU>
    if isempty(strfind(cell2mat(sheets), 'COM_Settings'))
        error('COM_Settings tab not in file')
    end
    [na1, na2, parameter] = xlsread(paramFile,'COM_Settings','',''); %#ok<ASGLU> % Import data from the settings file (imports the entire sheet)
catch ME %#ok<NASGU>
    warning('off','MATLAB:xlsread:Mode'); % suppress warning messages for reading the settings file from XLS
    [na1, na2, parameter] = xlsread(paramFile,'','','basic'); %#ok<ASGLU> % Import data from the settings file (imports the entire sheet)
end

% just need to define so we can pass
param.c=[.4e-12 .4e-12];
param.alen=[ 20 30 550  ];
param.az=[100 120 100];
% Default values are given for parameters when they are common to all clauses in 802.3bj and 803.2bm.

param.fb = xls_parameter(parameter, 'f_b')*1e9;
param.max_start_freq = xls_parameter(parameter, 'f_min')*1e9;
param.max_freq_step = xls_parameter(parameter, 'Delta_f')*1e9;
param.tx_ffe_c0_min = xls_parameter(parameter, 'c(0)', false);
param.tx_ffe_cm4_values = xls_parameter(parameter, 'c(-4)', true,0); % eval if string
param.tx_ffe_cm3_values = xls_parameter(parameter, 'c(-3)', true,0); % eval if string
param.tx_ffe_cm2_values = xls_parameter(parameter, 'c(-2)', true,0); % eval if string
param.tx_ffe_cm1_values = xls_parameter(parameter, 'c(-1)', true,0); % eval if string
param.tx_ffe_cp1_values = xls_parameter(parameter, 'c(1)', true,0); % eval if string
param.ndfe = xls_parameter(parameter, 'N_b');

param.g_DC_HP_values = xls_parameter(parameter, 'g_DC_HP', true,[]); % eval if string
param.f_HP = 1e9*xls_parameter(parameter, 'f_HP_PZ', true, []); % fp1 is in GHz
param.f_HP_Z = 1e9*xls_parameter(parameter, 'f_HP_Z', true, []); % fp1 is in GHz
param.f_HP_P = 1e9*xls_parameter(parameter, 'f_HP_P', true, []); % fp1 is in GHz

param.PHY_type= xls_parameter(parameter, 'PHY_type', false,'C2C'); % 
param.Min_VEO= xls_parameter(parameter, 'EH_min', true,[]); % for testing MTF and C2M
param.Max_VEO= xls_parameter(parameter, 'EH_max', true,[]); % for testing MTF and C2M

param.CTLE_type= xls_parameter(parameter, 'CTLE_type', false,'CL93'); % eval if string
if ~isempty(param.g_DC_HP_values) ; param.CTLE_type='CL120d';end % overrides CL93 if g_DC_HP_values is a spreadsheet entry
if ~isempty(param.f_HP_Z) ; param.CTLE_type='CL120e';end % overrides CL93 and CL120d if f_HP_Z is a spreadsheet entry
% always read in main ctle values. They would be interpreted different baseed
% on the clause they apply because of different CTF equations
param.ctle_gdc_values = xls_parameter(parameter, 'g_DC', true); % eval if string
param.CTLE_fp1 = 1e9*xls_parameter(parameter, 'f_p1', true, param.fb/4); % fp1 is in GHz
param.CTLE_fp2 = 1e9*xls_parameter(parameter, 'f_p2', true, param.fb); % fp2 is in GHz
param.CTLE_fz = 1e9*xls_parameter(parameter, 'f_z', true, param.fb/4); % fz is in GHz
% the contex of the poles an zeros are determined by the clause
switch param.CTLE_type
    case 'CL93'
        param.ctle_gdc_values = xls_parameter(parameter, 'g_DC', true); % eval if string
        param.CTLE_fp1 = 1e9*xls_parameter(parameter, 'f_p1', true, param.fb/4); % fp1 is in GHz
        param.CTLE_fp2 = 1e9*xls_parameter(parameter, 'f_p2', true, param.fb); % fp2 is in GHz
        param.CTLE_fz = 1e9*xls_parameter(parameter, 'f_z', true, param.fb/4); % fz is in GHz
    case 'CL120d'
        param.g_DC_HP_values = xls_parameter(parameter, 'g_DC_HP', true,[]); % eval if string
        param.f_HP = 1e9*xls_parameter(parameter, 'f_HP_PZ', true, []); % fp1 is in GHz
    case 'CL120e'
        % re adjust to get TD_CTLE to work with C:120e equation without
        % changing TD_CTLE code
        param.CTLE_fz =param.CTLE_fz ./ 10.^(param.ctle_gdc_values/20);
end

param.a_thru = xls_parameter(parameter, 'A_v', true);
param.a_fext = xls_parameter(parameter, 'A_fe', true);
param.a_next = xls_parameter(parameter, 'A_ne', true);
param.levels = xls_parameter(parameter, 'L');
param.specBER = xls_parameter(parameter, 'DER_0');
param.pass_threshold = xls_parameter(parameter, 'COM Pass threshold');
param.sigma_RJ = xls_parameter(parameter, 'sigma_RJ');
param.A_DD = xls_parameter(parameter, 'A_DD');
param.eta_0 = xls_parameter(parameter, 'eta_0');

param.SNDR = xls_parameter(parameter, 'SNR_TX', true);
param.R_LM = xls_parameter(parameter, 'R_LM');

param.samples_per_ui = xls_parameter(parameter, 'M', 32);
% This will keep bmax length 0 if Nb=0
param.bmax(1:param.ndfe) = xls_parameter(parameter, 'b_max(1)');
param.bmax(2:param.ndfe) = xls_parameter(parameter, 'b_max(2..N_b)', true);

% eval if string for all three - can use different for TX and RX
param.C_pkg_board = xls_parameter(parameter, 'C_p', true)*1e-9; % C_p in nF
param.C_diepad = xls_parameter(parameter, 'C_d', true)*1e-9; % C_d in nF
param.R_diepad = xls_parameter(parameter, 'R_d', true);
param.Z_t = xls_parameter(parameter, 'Z_t', true,50);
param.TR_TDR = xls_parameter(parameter, 'TR_TDR', true , 8e-3);


param.Z0 = xls_parameter(parameter, 'R_0', 50);
param.z_p_tx_cases = xls_parameter(parameter, 'z_p (TX)', true); % eval if string
param.z_p_next_cases = xls_parameter(parameter, 'z_p (NEXT)', true); % eval if string
param.z_p_fext_cases = xls_parameter(parameter, 'z_p (FEXT)', true); % eval if string
param.z_p_rx_cases = xls_parameter(parameter, 'z_p (RX)', true); % eval if string

% Table 93A-3 parameters
param.pkg_gamma0_a1_a2 = xls_parameter(parameter, 'package_tl_gamma0_a1_a2', true, [0 1.734e-3 1.455e-4]);
param.pkg_tau = xls_parameter(parameter, 'package_tl_tau', true, 6.141e-3);
param.pkg_Z_c = xls_parameter(parameter, 'package_Z_c', true, 78.2);

% Table 92-12 parameters
param.brd_gamma0_a1_a2 = xls_parameter(parameter, 'board_tl_gamma0_a1_a2', true, [0 4.114e-4 2.547e-4]); % eval string, default
param.brd_tau = xls_parameter(parameter, 'board_tl_tau', false, 6.191e-3);
param.brd_Z_c = xls_parameter(parameter, 'board_Z_c', true, 109.8);

param.z_bp_tx = xls_parameter(parameter, 'z_bp (TX)', false, 151);
param.z_bp_next = xls_parameter(parameter, 'z_bp (NEXT)', false, 72);
param.z_bp_fext = xls_parameter(parameter, 'z_bp (FEXT)', false, 72);
param.z_bp_rx = xls_parameter(parameter, 'z_bp (RX)', false, 151);

% Unofficial parameters
param.snpPortsOrder = xls_parameter(parameter, 'Port Order', true, [1 3 2 4]);
param.delta_IL=xls_parameter(parameter, 'delta_IL', false, 1);
% Deprecated parameters - affect only frequency domain analysis. better not change.
param.f_v = xls_parameter(parameter, 'f_v', false, 4);
param.f_f = xls_parameter(parameter, 'f_f', false, 4);
param.f_n = xls_parameter(parameter, 'f_n', false, 4);
param.f_r = xls_parameter(parameter, 'f_r', false, 4);
param.fb_BT_cutoff= xls_parameter(parameter, 'TDR_f_BT_3db', false, 20);
param.beta_x= xls_parameter(parameter, 'beta_x', false, 0);% for pkg adj 10e9
param.rho_x= xls_parameter(parameter, 'rho_x', false, .618); 
param.tfx= xls_parameter(parameter, 'fixture delay time', false, -1); ;
param.Grr_limit=xls_parameter(parameter, 'Grr_limit', false, 0);

% Operational control variables
%OP.include_pcb = xls_parameter(parameter, 'Include PCB (table 92-13)', false, 0);


OP.include_pcb = xls_parameter(parameter, 'Include PCB', false); % wrong table reference! use no default
OP.exit_if_deployed = xls_parameter(parameter, 'exit if deployed', false,0); % wrong table reference! use no default
OP.INCLUDE_CTLE = xls_parameter(parameter, 'INCLUDE_CTLE', false, 1);
OP.INCLUDE_FILTER = xls_parameter(parameter, 'INCLUDE_TX_RX_FILTER', false, 1);
OP.force_pdf_bin_size = xls_parameter(parameter, 'Force PDF bin size', false, 0);
OP.BinSize = xls_parameter(parameter, 'PDF bin size', false, 1e-5);
OP.DEBUG = xls_parameter(parameter, 'DIAGNOSTICS', false, false);
OP.DISPLAY_WINDOW = xls_parameter(parameter, 'DISPLAY_WINDOW', false, true);
OP.CSV_REPORT = xls_parameter(parameter, 'CSV_REPORT', false, true);
OP.SAVE_FIGURES=xls_parameter(parameter, 'SAVE_FIGURES', false, false);
OP.SAVE_FIGURE_to_CSV=xls_parameter(parameter, 'SAVE_FIGURE_to_CSV', false, false);
OP.GET_FD = xls_parameter(parameter, 'Display frequency domain', false, false);
OP.INC_PACKAGE = xls_parameter(parameter, 'INC_PACKAGE', false, true);
OP.IDEAL_RX_TERM = xls_parameter(parameter, 'IDEAL_RX_TERM', false, false);
OP.IDEAL_TX_TERM = xls_parameter(parameter, 'IDEAL_TX_TERM', false, false);
OP.EXTERNAL = xls_parameter(parameter, 'USE_EXTERNAL_PARAM',false, false);
OP.RESULT_DIR = regexprep(xls_parameter(parameter, 'RESULT_DIR'), '\\', filesep);
OP.RESULT_DIR=strrep(OP.RESULT_DIR,'{date}',date);
OP.BREAD_CRUMBS = xls_parameter(parameter, 'BREAD_CRUMBS',false, false);
OP.COM_CONTRIBUTION_CURVES = xls_parameter(parameter, 'COM_CONTRIBUTION',false,0);
OP.ENFORCE_CAUSALITY = xls_parameter(parameter, 'Enforce Causality', false, 0);
OP.EC_REL_TOL = xls_parameter(parameter, 'Enforce Causality REL_TOL', false, 1e-2);
OP.EC_DIFF_TOL = xls_parameter(parameter, 'Enforce Causality DIFF_TOL', false, 1e-3);
OP.EC_PULSE_TOL = xls_parameter(parameter, 'Enforce Causality pulse start tolerance', false, 0.01);
OP.pkg_len_select = xls_parameter(parameter, 'z_p select', true, 1);  % eval if string
OP.RX_CALIBRATION = xls_parameter(parameter, 'RX_CALIBRATION', false, false);
OP.sigma_bn_STEP = xls_parameter(parameter, 'Sigma BBN step', false, 5e-3);
OP.BBN_Q_factor = xls_parameter(parameter, 'BBN Q factor', false, 5);
OP.force_BBN_Q_factor = xls_parameter(parameter, 'Force BBN Q factor', false, false);
OP.transmitter_transition_time = xls_parameter(parameter, 'T_r', true , 8e-3);
OP.RL_norm_test=~xls_parameter(parameter, 'ERL_FOM', false, 0);

OP.T_r_meas_point = xls_parameter(parameter, 'T_r_meas_point', false, 0);
OP.T_r_filter_type= xls_parameter(parameter, 'T_r_filter_type', false, 0);
OP.FORCE_TR = xls_parameter(parameter, 'FORCE_TR', false, false);
% Control with OP.T_r_filter_type and OP.T_r_meas_point for backward
% compatibility
if OP.FORCE_TR
    OP.T_r_meas_point=0;
    OP.T_r_filter_type=1;
end
OP.TDR = xls_parameter(parameter, 'TDR', false, false);
OP.TDR_duration= xls_parameter(parameter, 'TDR_duration', false, 5);
OP.WC_PORTZ = xls_parameter(parameter, 'WC_PORTZ', false, false);
OP.T_k= xls_parameter(parameter, 'T_k', false, .6)*1e-9;
OP.ERL_ONLY = xls_parameter(parameter, 'ERL_ONLY', false,0); % compute ERL only 
OP.ERL=xls_parameter(parameter, 'ERL', false, false);
if OP.ERL , 
    OP.PTDR=1;
else
    OP.PTDR=0;
end % ERL needs to do a TDR
if OP.WC_PORTZ , OP.TDR=1;end % WC_PORTZ needs to do a TDR
OP.ADJ_4_LOAD = xls_parameter(parameter, 'ADJ_4_LOAD', false,0);
OP.BesselThompson=xls_parameter(parameter, 'BesselThompson', false, false);
OP.Butterworth=xls_parameter(parameter, 'TDR_Butterworth', false, true);
OP.inc_reflect_board=xls_parameter(parameter, 'inc_reflect_board', false,0);


OP.LIMIT_JITTER_CONTRIB_TO_DFE_SPAN = xls_parameter(parameter, 'LIMIT_JITTER_CONTRIB_TO_DFE_SPAN', false, false);
OP.impulse_response_truncation_threshold = xls_parameter(parameter, 'Impulse response truncatio threshold', false, 1e-3);
OP.interp_sparam_mag = xls_parameter(parameter, 'S-parameter magnitude extrapolation policy', false, 'linear_trend_to_DC');
OP.interp_sparam_phase = xls_parameter(parameter, 'S-parameter phase extrapolation policy', false, 'extrap_cubic_to_dc_linear_to_inf');
OP.PHY = xls_parameter(parameter, 'PHY', false, '');
OP.RUNTAG = xls_parameter(parameter, 'RUNTAG', false, '');
if isnan(OP.RUNTAG), OP.RUNTAG='';end
if isnumeric(OP.RUNTAG), OP.RUNTAG=num2str(OP.RUNTAG);end

% Parameters for error burst probability calculation. Not officially used.
OP.use_simple_EP_model = xls_parameter(parameter, 'Use simple error propagation model', false, false);
OP.nburst = xls_parameter(parameter, 'Max burst length calculated', false, 0);
OP.COM_EP_margin = xls_parameter(parameter, 'Error propagation COM margin', false, 0);


function p=xls_parameter(param_sheet, param_name, eval_if_string, default_value)
% helper function to read parameter values from XLS file. Uses names to find values.
if nargin<3, eval_if_string=0; end
[row, col]=find(strcmp(param_sheet, param_name));
if numel(row)*numel(col)==0
    if nargin<4
        missingParameter(param_name);
    else
        p = default_value;
    end
elseif numel(row)*numel(col)>1
    % if there are several occurrences, use the first, but warn
    warning('COM:XLS_parameter:MultipleOccurrence', ...
        '%d occurrences of "%s" found. Using the first', numel(row), param_name);
    p=param_sheet{row(1), col(1)+1};
else
    p=param_sheet{row, col+1};
end
if ischar(p) && eval_if_string
    p=eval(p);
end

function missingParameter (parameterName)
error( 'error:badParameterInformation', ...
    'The data for mandatory parameter %s is missing or incorrect' , parameterName);

function [sdd21p,param] = s21_pkg(sdd21, sdd11, sdd22, faxis, param, OP, channel_type, channel_number,chdatalocal)
% concatenates package reflections with s21,s11,and s22 with spec return loss (gammas)
% faxis is the frequency array
% sdd21, sdd11, sdd22 are the corresponding array of differential parameters

Z0=param.Z0;

% The following three parameters have possibly different valuesF for TX and
% RX (so can be 2-element vectors).
R_diepad = param.R_diepad;
C_diepad = param.C_diepad;
C_pkg_board = param.C_pkg_board;

% generate TX package according to channel type.
switch upper(channel_type)
    case 'THRU'
        pkg_len = param.Pkg_len_TX;
    case 'NEXT'
        pkg_len = param.Pkg_len_NEXT;
    case 'FEXT'
        pkg_len = param.Pkg_len_FEXT;
end
pkg_Z_Tx=param.pkg_Z_c(1);
pkg_Z_Rx=param.pkg_Z_c(2);
param.Tx_zc_sel=1;
param.Rx_zc_sel=2;
param.Tx_rd_sel=1; 
param.Rx_rd_sel=2;

% after all this param.Tx_seland param.Tx_sel are only need for VTF gamma
% we on need to set comming out if we have THRU channel

[ s11out, s12out, s21out, s22out ]= make_pkg(faxis, pkg_len, C_diepad(1), C_pkg_board(1),pkg_Z_Tx, param);

% RX package length is assumed to be the same for all channel types.
[ s11in, s12in, s21in, s22in ]= make_pkg(faxis, param.Pkg_len_RX, C_diepad(2), C_pkg_board(2), pkg_Z_Rx, param);

if OP.IDEAL_TX_TERM || (OP.RX_CALIBRATION == 1 && channel_number == 2) || OP.include_pcb == 2
    gamma_tx=0;
else
    gamma_tx=(R_diepad(param.Tx_rd_sel)-Z0)/(R_diepad(param.Tx_rd_sel)+Z0);% equation 93A-17
end
if OP.IDEAL_RX_TERM
    gamma_rx=0;
else
    gamma_rx=(R_diepad(param.Rx_rd_sel)-Z0)/(R_diepad(param.Rx_rd_sel)+Z0);% equation 93A-17
end

if OP.INC_PACKAGE==0
    sdd21p= sdd21;
else
    if OP.RX_CALIBRATION == 1 && channel_number == 2 
        %   for calibration do not include the transmitter package
        [s11out_rx, UNUSED_OUTPUT, s21out_rx, s22out_rx ] = combines4p( sdd11, sdd21, sdd21, sdd22, s22out, s12out, s21out, s11out ); %#ok<ASGLU> % s22 is ball side of package
        %% Equation 93A-18
        sdd21p= s21out_rx.*(1-gamma_tx).*(1+gamma_rx)./(1.- s11out_rx.*gamma_tx - s22out_rx.*gamma_rx  -s21out_rx.^2.*gamma_tx.*gamma_rx +s11out_rx.*s22out_rx.*gamma_tx.*gamma_rx);
    else
        %% Equations 93A-4 to 93A-7
        if ~OP.IDEAL_TX_TERM 
            [sdd11, UNUSED_OUTPUT, sdd21, sdd22] = combines4p( s11out, s12out, s21out, s22out, sdd11, sdd21, sdd21, sdd22 ); %#ok<ASGLU>
        end
        if OP.IDEAL_TX_TERM ||  OP.T_r_filter_type == 1 
         % for RITT testing with good termination as in some instruments
         % and tx filter when required
            if OP.T_r_filter_type==0
                H_t = exp(-(pi*faxis/1e9*OP.transmitter_transition_time/1.6832).^2); %% Equation 93A-46 %%
            else
                tr=OP.transmitter_transition_time;
                f9=faxis/1e9;
                if OP.T_r_meas_point == 1
                    k=1.9466+7.12*sqrt(1-6.51e-3/tr);
                    H_t=105./(f9.^4*(k*tr)^4 - f9.^3*(k*tr)^3*10i - 45*f9.^2*(k*tr)^2 + f9*(k*tr)*105i + 105);
                else
                  H_t = exp( -2*(pi*f9*tr/1.6832).^2 ).*exp(-1j*2*pi*f9*tr*3);
                 end
                                        
            end
            sdd21 = sdd21 .* H_t;
        end
        if ~OP.IDEAL_RX_TERM
            [sdd11, UNUSED_OUTPUT, sdd21, sdd22] = combines4p( sdd11, sdd21, sdd21, sdd22, s22in, s12in, s21in, s11in ); %#ok<ASGLU> % s22 is ball side of package
        end
        %% Equation 93A-18
        sdd21p= sdd21.*(1-gamma_tx).*(1+gamma_rx)./(1.- sdd11.*gamma_tx - sdd22.*gamma_rx  -sdd21.^2.*gamma_tx.*gamma_rx +sdd11.*sdd22.*gamma_tx.*gamma_rx);
    end
end


function [ s11out, s12out, s21out, s22out ] = make_pkg(f, pkg_len, cpad, cball, pkg_z, param)
f(f<eps)=eps;
tau = param.pkg_tau; gamma0_a1_a2=param.pkg_gamma0_a1_a2; Lenscale=pkg_len; zref=param.Z0;
%% Equation 93A-8
s11pad= -1i*2*pi.*f*cpad*zref./(2+1i*2*pi.*f*cpad*zref);
s21pad= 2./(2+1i*2*pi.*f*cpad*zref);

[ S11, S12, S21, S22 ] = synth_tline(f, pkg_z, zref, gamma0_a1_a2, tau, Lenscale); %#ok<NASGU,ASGLU>
[ s11out1, s12out1, s21out1, s22out1 ]= ...
    combines4p(  s11pad, s21pad, s21pad, s11pad, S11, S21, S21, S11 ); % first part of equation 93A-15

%% Equation 93A-8
s11ball= -1i*2*pi.*f*cball*zref./(2+1i*2*pi.*f*cball*zref);
s21ball= 2./(2+1i*2*pi.*f*cball*zref);
[ s11out, s12out, s21out, s22out ]= ...
    combines4p( s11out1, s12out1, s21out1, s22out1, s11ball, s21ball, s21ball, s11ball );% second part of equation 93A-15

function [ s11out, s12out, s21out, s22out ] = add_brd(chdata, param, OP)
% Used in Clause 92 for adding board trace between TP0 and TP2
switch chdata.type
    case 'THRU'
        z_bp_tx = param.z_bp_tx;
    case 'NEXT'
        z_bp_tx = param.z_bp_next;
    case 'FEXT'
        z_bp_tx = param.z_bp_fext;
end

[ s11tx, s12tx, s21tx, s22tx ] = synth_tline(chdata.faxis, param.brd_Z_c(1), param.Z0, param.brd_gamma0_a1_a2, param.brd_tau, z_bp_tx);
[ s11rx, s12rx, s21rx, s22rx ] = synth_tline(chdata.faxis, param.brd_Z_c(2), param.Z0, param.brd_gamma0_a1_a2, param.brd_tau, param.z_bp_rx);

switch OP.include_pcb
    case 1
        [ s11out1, s12out1, s21out1, s22out1 ]=combines4p(  s11tx, s12tx, s21tx, s22tx, chdata.sdd11_raw, chdata.sdd21_raw, chdata.sdd12_raw, chdata.sdd22_raw );
        [ s11out, s12out, s21out, s22out ]=combines4p(  s11out1, s12out1, s21out1, s22out1,  s11rx, s12rx, s21rx, s22rx);
    case 2
        [ s11out, s12out, s21out, s22out ]=combines4p( chdata.sdd11_raw, chdata.sdd21_raw, chdata.sdd12_raw, chdata.sdd22_raw, s11rx, s12rx, s21rx, s22rx);    
end

function [ s11out, s12out, s21out, s22out ] = add_brdexp(chdata, param, OP)
% Used in Clause 92 for adding board trace between TP0 and TP2
switch chdata.type
    case 'THRU'
        z_bp_tx = param.z_bp_tx;
    case 'NEXT'
        z_bp_tx = param.z_bp_next;
    case 'FEXT'
        z_bp_tx = param.z_bp_fext;
end
zref=50;
c1=1.8e-13;
c2=1.8e-13;
f=chdata.faxis;
f(f<eps)=eps;
s11pad1= -1i*2*pi.*f*c1*zref./(2+1i*2*pi.*f*c1*zref);
s21pad1= 2./(2+1i*2*pi.*f*c1*zref);
s11pad2= -1i*2*pi.*f*c2*zref./(2+1i*2*pi.*f*c2*zref);
s21pad2= 2./(2+1i*2*pi.*f*c2*zref);
param.Zp_adder=100; % mm
param.brd_Z_c(3)=180; %ohms
% add lead in
[ s11txa, s12txa, s21txa, s22txa ] = synth_tline(chdata.faxis, param.brd_Z_c(3), param.Z0, param.brd_gamma0_a1_a2, param.brd_tau, param.Zp_adder);
% create line at tx end
[ s11tx, s12tx, s21tx, s22tx ] = synth_tline(chdata.faxis, param.brd_Z_c(1), param.Z0, param.brd_gamma0_a1_a2, param.brd_tau, z_bp_tx);
[s11tx, s12tx, s21tx, s22tx  ]= ...
    combines4p(  s11xa, s21xa, s21xa, s11xa, s11tx, s12tx, s21tx, s22tx   ); 


% surround with tx voar caps
[s11tx, s12tx, s21tx, s22tx  ]= ...
    combines4p(  s11pad1, s21pad1, s21pad1, s11pad1, s11tx, s12tx, s21tx, s22tx   ); 
[s11tx, s12tx, s21tx, s22tx  ]= ...
    combines4p(   s11tx, s12tx, s21tx, s22tx,s11pad2, s21pad2, s21pad2, s11pad2  );


[ s11rx, s12rx, s21rx, s22rx ] = synth_tline(chdata.faxis, param.brd_Z_c(2), param.Z0, param.brd_gamma0_a1_a2, param.brd_tau, param.z_bp_rx);

switch OP.include_pcb
    case 1
        [ s11out1, s12out1, s21out1, s22out1 ]=combines4p(  s11tx, s12tx, s21tx, s22tx, chdata.sdd11_raw, chdata.sdd21_raw, chdata.sdd12_raw, chdata.sdd22_raw );
        [ s11out, s12out, s21out, s22out ]=combines4p(  s11out1, s12out1, s21out1, s22out1,  s11rx, s12rx, s21rx, s22rx);
    case 2
        [ s11out, s12out, s21out, s22out ]=combines4p( chdata.sdd11_raw, chdata.sdd21_raw, chdata.sdd12_raw, chdata.sdd22_raw, s11rx, s12rx, s21rx, s22rx);    
end

function [s11, s12, s21, s22] = synth_tline(f, Z_c, Z_0, gamma_coeff, tau, d)
f_GHz=f/1e9;
%% Equation 93A-10 %%
gamma_1 = gamma_coeff(2)*(1+1i);
%% Equation 93A-11 %%
gamma_2 = gamma_coeff(3)*(1-2i/pi*log(f_GHz)) + 2i*pi*tau;
%% Equation 93A-9 %%
gamma = gamma_coeff(1)+gamma_1.*sqrt(f_GHz)+gamma_2.*f_GHz;
gamma(f_GHz==0) = gamma_coeff(1);

%% Equation 93A-12 %%
rho_rl=(Z_c-2*Z_0)/(Z_c+2*Z_0);

exp_gamma_d = exp(-d*gamma);
%% Equations 93A-13 and 93A-14 %%
s11 = rho_rl*(1-exp_gamma_d.^2)./(1-rho_rl^2*exp_gamma_d.^2);
s21 = (1-rho_rl^2)*exp_gamma_d./(1-rho_rl^2*exp_gamma_d.^2);
s12 = s21;
s22 = s11;

function [ s11out, s12out, s21out, s22out ] = combines4p( s11in1, s12in1, s21in1, s22in1, s11in2, s12in2, s21in2, s22in2)

s1=zeros(2,2,length(s11in1)); s2=s1; t3=s1;
for i=1:length(s11in1)
    s1(:,:,i)=[s11in1(i) s12in1(i); s21in1(i) s22in1(i) ];
    s2(:,:,i)=[s11in2(i) s12in2(i); s21in2(i) s22in2(i) ];
end
t1=stot(s1);
t2=stot(s2);
for i=1:length(s11in1)
    t3(:,:,i)=t1(:,:,i)*t2(:,:,i);
end
s3=ttos(t3);
s11out=s3(1,1,:);
s11out=transpose(s11out(:));
s12out=s3(1,2,:);
s12out=transpose(s12out(:));
s21out=s3(2,1,:);
s21out=transpose(s21out(:));
s22out=s3(2,2,:);
s22out=transpose(s22out(:));

function t_params = stot(s_params)
% p 67 R. Mavaddat. (1996). Network scattering parameter. Singapore: World Scientific.
% ISBN 978-981-02-2305-2. http://books.google.com/?id=287g2NkRYxUC&lpg=PA65&dq=T-parameters+&pg=PA67.
[s11, s12, s21, s22] = deal(s_params(1,1,:), s_params(1,2,:), s_params(2,1,:), s_params(2,2,:));
delta = (s11.*s22-s12.*s21);
s21(s21==0)=eps;
t_params = [1./s21, -s22./s21; s11./s21, -delta./s21];

function s_params = ttos(t_params)
% p 67 R. Mavaddat. (1996). Network scattering parameter. Singapore: World Scientific.
% ISBN 978-981-02-2305-2. http://books.google.com/?id=287g2NkRYxUC&lpg=PA65&dq=T-parameters+&pg=PA67.
[t11, t12, t21, t22] = deal(t_params(1,1,:), t_params(1,2,:), t_params(2,1,:), t_params(2,2,:));
delta = t11.*t22-t21.*t12;
t11(t11==0)=eps;
s_params = [t21./t11, delta./t11; 1./t11, -t12./t11];

function [ sigma_NE ] = get_sigma_noise( H_ctf, param, chdata, sigma_bn )
% the FEXT channel for calibration basically a DC connection unlike normal
% FEXT channels which are nearly open at DC channels
H_r = 1./polyval([1 2.613126 3.414214 2.613126 1], 1i*chdata(2).faxis/(param.f_r*param.fb));
idxfbby2=find( chdata(2).faxis(:) >= param.fb/2, 1);
if size(chdata,2) >= 2
    Hnoise_channel=chdata(2).sdd21;% rx package is already included tx is not
else
    Hnoise_channel=1;
end
H_np=Hnoise_channel.*H_ctf.*H_r;
%% Equation 93A-48 %%
sigma_NE = sigma_bn*sqrt(mean(abs(H_np(1:idxfbby2).^2)));

function sigma_XT = get_xtlk_noise( upsampled_txffe, type, param, chdata )
% calculate crosstalk sigma at worst phase, per equation 93A-33, for
% channels of chosen type (FEXT and NEXT are treated separately).

sigma_XT_sqd=0;
for k=1:param.num_s4p_files
    if isequal(chdata(k).type, type)
        sigma_i_sqd=0;
        effective_channel = filter(upsampled_txffe, 1, chdata(k).ctle_imp_response);
        N=round(length(effective_channel)/param.samples_per_ui)-3;
        sbr=filter(ones(param.samples_per_ui, 1), 1, effective_channel);
        for m=1:param.samples_per_ui
            h=sbr(param.samples_per_ui*(1:N)+m);
            sigma_i_sqd=max(param.sigma_X^2*sum(h.^2),sigma_i_sqd);
        end
        sigma_XT_sqd = sigma_XT_sqd + sigma_i_sqd;
    end
end
sigma_XT=sqrt(sigma_XT_sqd);

function [chdata, param] = get_s4p_files(param, OP, num_fext, num_next, file_list)
% filename parsing and acquisition
%------------------------------------------------------------------
%----------put files names into chdata structure ---------
% The thru file has the index of 1
% crosstalk file are indexed from 2
% nxi is incremented each time a file is read in  so that nxi will end
filepath=[]; % path name for file
nxi=0; % file index
% get the THRU file
if size(file_list,2) ~= 0
    file_list(1)=strrep(file_list(1),'\', filesep); % OS file convention conversion
    [filepath, basename, fileext]=fileparts(file_list{1});
    
else
    if OP.RX_CALIBRATION == 1 || OP.FIXTURE_CALIBRATION
        h = msgbox('enter test channel s-parameter file'); set(h,'Color',[1 .85 0]);
        movegui(h,'northeast')
    end
    if  OP.ERL == 2;
        dir=fullfile(filepath, '*.s2p');
        [basename,filepath]=uigetfile(dir,'input RL measurement .s2p ');
        if filepath == 0
            error('No RL measurement file')
        end
    else
        dir=fullfile(filepath, '*.s4p');
        [basename,filepath]=uigetfile(dir,'input thru channel response .s4p ');
        if filepath == 0
            error('No Thru file')
        end
    end
    [UNUSED_OUTPUT, basename, fileext]=fileparts(basename); %#ok<ASGLU>
end
nxi=nxi+1;
chdata(nxi).filename = fullfile(filepath,  [basename fileext]);
chdata(nxi).ext = fileext;
[UNUSED_OUTPUT, dirname]=fileparts(filepath); %#ok<ASGLU>
%    chdata(nxi).base=[pth(max(strfind(pth,filesep))+1:end) '--' basename]; % add 1 directory back to basename
chdata(nxi).base=[OP.RUNTAG ' ' dirname  '--' basename]; % add 1 directory back to basename
% chdata(nxi).A=param.A_thru; % pam encoding amplitude reduction is do in reporting but is comprehending in crosstalk PDF
chdata(nxi).type='THRU';
chdata(nxi).ftr=param.fb*param.f_v;
param.base=chdata(nxi).base; %for print out function that don't pass chdata

% now get FEXT file names into chdata structure
kxi=nxi;
for nxi=kxi+1:num_fext+kxi
    lastfilepath=filepath;
    if size(file_list,2) ~= 0
        [filepath, basename, fileext]=fileparts(file_list{nxi});
    else
        if OP.RX_CALIBRATION == 1 || OP.FIXTURE_CALIBRATION
            h = msgbox('enter noise channel s-parameter file'); set(h,'Color',[1 .85 0]);
            movegui(h,'northeast')
        end
        if  param.tfx == -1 && OP.ERL_ONLY && OP.ERL == 2;
            dir=fullfile(filepath, '*.s4p');
                [basename,filepath]=uigetfile(dir,'input noise channel response .s4p');
            if filepath==0
                error('Not enough NEXT files')
            end
        else
            dir=fullfile(filepath, '*.s4p');
            if OP.RX_CALIBRATION == 1 || OP.FIXTURE_CALIBRATION
                [basename,filepath]=uigetfile(dir,'input noise channel response .s4p');
            else
                [basename,filepath]=uigetfile(dir,['input fext channel response .s4p #', num2str(nxi-kxi)]);
            end
            if filepath==0
                error('Not enough NEXT files')
            end
        end
        [UNUSED_OUTPUT, basename, fileext]=fileparts(basename); %#ok<ASGLU>
    end
    if isempty( filepath), filepath=lastfilepath; end
    chdata(nxi).filename = fullfile(filepath,  [basename fileext]);
    chdata(nxi).ext = fileext;
    [UNUSED_OUTPUT, dirname]=fileparts(filepath); %#ok<ASGLU>
    chdata(nxi).base=[OP.RUNTAG ' ' dirname '--' basename ];
    %     chdata(nxi).A=param.a_fext;
    chdata(nxi).ftr=param.fb*param.f_f;
    chdata(nxi).type='FEXT';
end
% now get NEXT file names into chdata structure
kxi=num_fext+kxi;
for nxi=kxi+1:num_next+kxi
    lastfilepath=filepath;
    if size(file_list,2) ~= 0
        [filepath, basename, fileext]=fileparts(file_list{nxi});
    else
        dir=fullfile(filepath, '*.s4p');
        [basename,filepath]=uigetfile(dir,['input next channel response .s4p ', num2str(nxi-kxi)]);
        if filepath==0
            error('Not enough NEXT files')
        end
        [UNUSED_OUTPUT, basename, fileext]=fileparts(basename); %#ok<ASGLU>
    end
    if isempty( filepath), filepath=lastfilepath; end
    chdata(nxi).filename = fullfile(filepath,  [basename fileext]);
    chdata(nxi).ext = fileext;
    [UNUSED_OUTPUT, dirname]=fileparts(filepath); %#ok<ASGLU>
    chdata(nxi).base=[OP.RUNTAG ' ' dirname '--' basename ];
    %     chdata(nxi).A=param.A_next;
    chdata(nxi).ftr=param.fb*param.f_n;
    chdata(nxi).type='NEXT';
end

function [chdata, param, SDDch ] = read_s4p_files(param, OP, chdata)
%% extract s-parameter and convert to differential mode
% extract s-parameter data from files and apply tx and rx filters as well as package filters
num_files=length(chdata);
if ~OP.DISPLAY_WINDOW, fprintf('reading file '); end
for i=1:num_files
    if OP.DISPLAY_WINDOW; hwaitbar=waitbar(0);end
    progress = i/num_files;
    if OP.DISPLAY_WINDOW
        waitbar(progress, hwaitbar, ['Processing ' chdata(i).base]); figure(hwaitbar); drawnow;
    else
        fprintf('%i ',i);
    end
    
    % Skip reading file if it was already read (multiple test cases)
    if (~isfield(chdata(i), 'faxis')) || isempty(chdata(i).faxis)
        param.FLAG.S2P=0;
        switch chdata(i).ext
            case '.s2p' % for differential return loss
                [Sch,SDDch] = read_p2_s2params(chdata(i).filename,  0, 0, param.snpPortsOrder, OP);
                chdata(i).fmaxi = length(Sch.freq);
                chdata(i).faxis = Sch.freq;
                chdata(i).sdd11_raw = transpose(SDDch(1:chdata(i).fmaxi,1,1));
                param.FLAG.S2P=1;
            case '.s4p'
                if ~isfield(param,'snpPortsOrder') || isempty(param.snpPortsOrder)
                    param.snpPortsOrder = [1 3 2 4]; % default order normally used.
                elseif length(param.snpPortsOrder) ~= 4
                    error( 'warning:sNpFilePortMismatch', ...
                        '\n\t The number of ports defined (%G) does not match the sNp file type (%s)', ...
                        length(param.snpPortsOrder), ...
                        chdata(i).ext ...
                        );
                end
                % read function returns differnetial mode parameters
                [Sch, SDDch] = read_p4_s4params(chdata(i).filename,  0, 0, param.snpPortsOrder, OP);
                chdata(i).fmaxi = length(Sch.freq);
                
                if Sch.freq(chdata(i).fmaxi) < param.fb
                    warning('COM:read_s4p:MaxFreqTooLow', ...
                        'In %s: the maximum frequency provided, %g, is less than the signaling rate: %g', ...
                        chdata(i).filename, Sch.freq(end), param.fb);
                end
                if Sch.freq(1) > param.max_start_freq
                    warning('COM:read_s4p:StartFreqTooHigh', ...
                        'In %s: minimum frequency, %.2g GHz, is larger than the recommended %.2g GHz', ...
                        chdata(i).filename, Sch.freq(1)/1e9, param.max_start_freq/1e9);
                end
                freqstep=diff(Sch.freq);
                % ignore frequency differences up to 1 Hz - possible numerical artifacts
                if max(freqstep)-min(freqstep) > 1
                    warning('COM:read_s4p:NonUniformFreqSpacing', 'In %s: non-uniform frequency steps: min=%.3g GHz, max=%.3g GHz', ...
                        chdata(i).filename, min(freqstep)/1e9, max(freqstep)/1e9);
                end
                if max(freqstep) - param.max_freq_step > 1
                    warning('COM:read_s4p:FreqStepTooHigh', 'In %s: frequency step, %.2g GHz, is larger than the recommended %.2g GHz', ...
                        chdata(i).filename, max(freqstep)/1e9, param.max_freq_step/1e9);
                end
                
                chdata(i).faxis = Sch.freq;
                chdata(i).sdd12_raw = transpose(SDDch(1:chdata(i).fmaxi,1,2));
                chdata(i).sdd21_raw = transpose(SDDch(1:chdata(i).fmaxi,2,1));
                chdata(i).sdd22_raw = transpose(SDDch(1:chdata(i).fmaxi,2,2));
                chdata(i).sdd11_raw = transpose(SDDch(1:chdata(i).fmaxi,1,1));
            otherwise
        end
        %
        % differential response is sdd21 <------------------------------
        % differential return losses are sdd11 and sdd22
        % save original
        switch chdata(i).ext
            case '.s2p' % for differential return loss
                chdata(i).sdd11_orig=chdata(i).sdd11_raw;
                chdata(i).sdd11=chdata(i).sdd11_raw;
            otherwise
                chdata(i).sdd11_orig=chdata(i).sdd11_raw;
                chdata(i).sdd22_orig=chdata(i).sdd22_raw;
                chdata(i).sdd12_orig=chdata(i).sdd12_raw;
                chdata(i).sdd21_orig=chdata(i).sdd21_raw;
                if OP.include_pcb
                    % add boards to sdd
                    [chdata(i).sdd11_raw, chdata(i).sdd12_raw, chdata(i).sdd21_raw, chdata(i).sdd22_raw] = add_brd(chdata(i), param, OP);
                    
                end
                chdata(i).sdd11=chdata(i).sdd11_raw;
                chdata(i).sdd22=chdata(i).sdd22_raw;
        end
    else
        SDDch(:,1,2)=chdata(i).sdd12_raw;
        SDDch(:,2,1)=chdata(i).sdd21_raw;
        SDDch(:,1,1)=chdata(i).sdd11_raw;
        SDDch(:,2,2)=chdata(i).sdd22_raw;
    end
    if ~param.FLAG.S2P
        if  OP.INC_PACKAGE ~= 0
            % the new sdd
            % changes to s21_pkg for calibration is to not include Tx package in the sp21p
            [chdata(i).sdd21p,param]= s21_pkg(chdata(i).sdd21_raw, chdata(i).sdd11_raw, chdata(i).sdd22_raw, chdata(i).faxis, param, OP, chdata(i).type, i, chdata(i));
            chdata(i).sdd21=chdata(i).sdd21p;
        else
            chdata(i).sdd21=chdata(i).sdd21_raw;
        end
        chdata(i).sdd21f=chdata(i).sdd21_raw; % used for FD analysis i.e. not filtered
    end
end
if ~OP.DISPLAY_WINDOW, fprintf('\n'); end

function plot_bathtub_curves(hax, max_signal, sci_pdf, cci_pdf, isi_and_xtalk_pdf, noise_pdf, combined_interference_and_noise_pdf, bin_size)
cursors = d_cpdf(bin_size,max_signal*[-1 1], [1 1]/2);
signal_and_isi_pdf = conv_fct(cursors, sci_pdf);
signal_and_xtalk_pdf = conv_fct(cursors, cci_pdf);
signal_and_channel_noise_pdf = conv_fct(cursors, isi_and_xtalk_pdf);
signal_and_system_noise_pdf = conv_fct(cursors, noise_pdf);
signal_and_total_noise_pdf = conv_fct(cursors, combined_interference_and_noise_pdf);
%% Added by Bill Kirkland, June 14, 2017
cursors_l = cursors; cursors_l.y(cursors_l.x>0) = 0;
cursors_r = cursors; cursors_r.y(cursors_r.x<0) = 0;
signal_and_total_noise_pdf_l = conv_fct(cursors_l, combined_interference_and_noise_pdf);
signal_and_total_noise_pdf_r = conv_fct(cursors_r, combined_interference_and_noise_pdf);

semilogy(signal_and_isi_pdf.x, abs(cumsum(signal_and_isi_pdf.y)-0.5) ,'r','Disp','ISI', 'parent', hax)
hold on
semilogy(signal_and_xtalk_pdf.x, abs(cumsum(signal_and_xtalk_pdf.y)-0.5) ,'b','Disp','Xtalk', 'parent', hax)
semilogy(signal_and_channel_noise_pdf.x, abs(cumsum(signal_and_channel_noise_pdf.y)-0.5) ,'c','Disp','ISI+Xtalk', 'parent', hax)
semilogy(signal_and_system_noise_pdf.x, abs(cumsum(signal_and_system_noise_pdf.y)-0.5) ,'m','Disp','Jitter, TX and system noise', 'parent', hax)
%% Added by Bill Kirkland, June 14, 2017
% modification allows bathtub curves to cross over and hence one can
% directly read the noise component.
%semilogy(signal_and_total_noise_pdf.x, abs(cumsum(signal_and_total_noise_pdf.y)-0.5) ,'k','Disp','total noise PDF', 'parent', hax)
vbt_l = abs(0.5-cumsum(signal_and_total_noise_pdf_l.y));
vbt_r = fliplr(0.5-(cumsum(fliplr(signal_and_total_noise_pdf_r.y))));
semilogy(signal_and_total_noise_pdf_l.x, vbt_l ,'k','Disp','total noise PDF left', 'parent', hax)
semilogy(signal_and_total_noise_pdf_r.x, vbt_r ,'k','Disp','total noise PDF right', 'parent', hax)

hc=semilogy(max_signal*[-1 -1 1 1], [0.5 1e-20 1e-20 0.5], '--ok');
set(get(get(hc,'Annotation'),'LegendInformation'), 'IconDisplayStyle','off');

ylabel(hax, 'Probability')
xlabel(hax, 'volts')
legend(hax, 'show')


function plot_pie_com(hax, max_signal, sci_pdf, cci_pdf, isi_and_xtalk_pdf, noise_pdf, combined_interference_and_noise_pdf, bin_size,param)
BER=param.specBER;
delta_dB=param.delta_IL;

iex.combined_interference_and_noise_pdf=find(abs(cumsum(combined_interference_and_noise_pdf.y))>= BER, 1, 'first');
iex.noise_pdf=find(abs(cumsum(noise_pdf.y))>= BER, 1, 'first');
iex.cci_pdf=find(abs(cumsum(cci_pdf.y))>= BER, 1, 'first');
iex.sci_pdf=find(abs(cumsum(sci_pdf.y))>= BER, 1, 'first');


maxn(1)=abs(sci_pdf.x(iex.sci_pdf));
maxn(2)=abs(noise_pdf.x(iex.noise_pdf));
maxn(3)=abs(cci_pdf.x(iex.cci_pdf));
maxn_tot=abs(combined_interference_and_noise_pdf.x(iex.combined_interference_and_noise_pdf));

COM=20*log10(max_signal/maxn_tot);
COM_per_noise(1:3)=COM*(maxn(1:3).^2/sum(maxn.^2));
COM_ISI=sprintf(  '%.2g%%',100*COM_per_noise(1)/sum(COM_per_noise));
COM_SYS=sprintf(  '%.2g%%',100*COM_per_noise(2)/sum(COM_per_noise));
COM_xtalk=sprintf('%.2g%%',100*COM_per_noise(3)/sum(COM_per_noise));

pfctr=exp(-0.09054*delta_dB);% less loss
mfctr=exp(0.09054*delta_dB); % more loss

%less loss
plus_maxn(1)=abs(sci_pdf.x(iex.sci_pdf))*pfctr;
plus_maxn(2)=abs(noise_pdf.x(iex.noise_pdf));
plus_maxn(3)=abs(cci_pdf.x(iex.cci_pdf))*pfctr;
% plus_maxn_tot=(maxn(1)*pfctr+maxn(2)+maxn(3)*pfctr)*maxn_tot/sum(plus_maxn);
plus_maxn_tot=norm(plus_maxn);

minus_maxn(1)=abs(sci_pdf.x(iex.sci_pdf))*mfctr;
minus_maxn(2)=abs(noise_pdf.x(iex.noise_pdf));
minus_maxn(3)=abs(cci_pdf.x(iex.cci_pdf))*mfctr;
% minus_maxn_tot=(maxn(1)*mfctr+maxn(2)+maxn(3)*mfctr)*maxn_tot/sum(minus_maxn);
minus_maxn_tot=norm(minus_maxn);

% more loss
COMp=20*log10(max_signal*pfctr/maxn_tot);
COMp_per_noise(1:3)=COMp*(plus_maxn(1:3).^2/plus_maxn_tot^2);
COMp_ISI=sprintf(  '%.2gdB',COMp_per_noise(1));
COMp_SYS=sprintf(  '%.2gdB',COMp_per_noise(2));
COMp_xtalk=sprintf('%.2gdB',COMp_per_noise(3));
% less loss
COMm=20*log10(max_signal*mfctr/maxn_tot);
COMm_per_noise(1:3)=COMm*(minus_maxn(1:3).^2/minus_maxn_tot^2);
COMm_ISI=sprintf(  '%.2gdB',COMm_per_noise(1));
COMm_SYS=sprintf(  '%.2gdB',COMm_per_noise(2));
COMm_xtalk=sprintf('%.2gdB',COMm_per_noise(3));

xax=[[delta_dB delta_dB delta_dB];[ 0 0 0 ]; [ -delta_dB -delta_dB -delta_dB]];


if(COM<0)
    return
end

labels= { ['ISI ' COM_ISI] ['System noise/jitter ' COM_SYS] [' Crosstalk ' COM_xtalk]};

% pie(COM_per_noise,labels)
% legend(labels,{'ISI COM','System noise/jitter COM','Crosstalk COM'});
% legend('show','Location','bestoutside')
nullbar= [ 0 0 0 ];
bar(xax,[nullbar ;COM_per_noise; nullbar],'stacked')
% bar(xax,[COMp_per_noise;COM_per_noise;COMm_per_noise],'stacked')
hold on
bar(xax,[[COMp 0 0 ];nullbar;nullbar],'stacked','w','barwidth',.5)
bar(xax,[nullbar ;nullbar; [0 0 COMm]],'stacked','w','barwidth',.5)

plot([-delta_dB,0, delta_dB], [param.pass_threshold ,param.pass_threshold ,param.pass_threshold ],'r')
% ax=gca;
% ax.XTickLabels = {'decreasing loss','-','increasing loss'};
set(gca,'XTickLabel', {'decreasing loss','-','increasing loss'})
grid on
legend(labels,'Location','north')
xlabel(['if loss is changed by ' sprintf('%.2gdB',delta_dB) ])
ylabel('COM (dB)')
hold off




function recolor_plots(ax)

if ~verLessThan('matlab', '8.4.0')
    return
end
colors='brgcmk';
ch=flipud(get(ax, 'children'));

for k=1:length(ch)
    set(ch(k), 'Color', colors(mod(k-1, length(colors))+1));
    set(ch(k), 'LineWidth', 2*floor((k-1)/length(colors))+1);
end
legend (ax, 'off');
warning('off', 'MATLAB:legend:PlotEmpty');
set(legend (ax, 'show'), 'interp', 'none');

function csv_string = str2csv(c)
% convert a cell array of strings to a csv string
cell_tmp = cell(2, length(c));
cell_tmp(1,:)=c;
cell_tmp(2,:) = {','};
cell_tmp{2,end} = '';
csv_string=strcat(cell_tmp{:});

function [ h ] = savefigs( param, OP )

%% find the figures
hw = waitbar(0,'Saving figures...');
h = findobj(0, 'Type', 'figure');
for ii=1:length(h)
    
    figname= get(h(ii), 'Name'); % use the figure name as file name
    if isempty(strfind(figname,param.base))
        figname = [figname ' ' OP.RUNTAG ' ' param.base ]; %#ok<AGROW>
    end
    if verLessThan('matlab', '8.4.0')
        figname = ['f_' num2str(h(ii)) '_' figname]; %#ok<AGROW>
    else
        figname = ['f_' num2str(h(ii).Number) '_' figname]; %#ok<AGROW>
    end
    figname = strrep(figname,':','-');
    figname = strrep(figname,' ','_');
    if OP.SAVE_FIGURES==1
        saveas(h(ii), fullfile(OP.RESULT_DIR, [figname '.fig']));
    end
    %% get x y data
    if OP.SAVE_FIGURE_to_CSV==1
        h_L = findobj(h(ii),'Type','line'); % find handles to all the lines
        M=[]; %ncol=1;
        for nk=1:length(h_L)
            % get x and data for a line.
            x_data=get(h_L(nk),'xdata')';
            y_data=get(h_L(nk),'ydata')';
            % .........>> need to get data in the line structure (legend or label) for headers
            M=[M; x_data; y_data]; %#ok<AGROW>
        end
        csvwrite([OP.RESULT_DIR figname '.csv'],M);
        %      clear M y x header h_L
    end
    waitbar(ii/length(h),hw)

end

close(hw)

%%
function TDR_results = get_TDR(sdd, OP, param,ZT)
% sdd is differential s-parameters structure (2 port assumed)
% input parameter structure for s parameters sdd--> sdd.Impedance, sdd.Frequencies, sdd.Parameters, sdd.NumPorts
% TDR_results.delay             pre t=0 delay for TDR... help with time domain responce quaility
% TDR_results.tdr               the TDR responce (ohms vs  TDR_results.t
% TDR_results.t                 starting at t=0
% TDR_results.tx_filter         transmitter filter vs TDR_results.f
% TDR_results.Rx_filter         receiver filter vs TDR_results.f
% TDR_results.f                 frequency for filter and s parameters
% TDR_results.ptdr_RL           reflection waveform from the pulse
% TDR_results.WC_ptdr_samples_t worst case time sample of the reflection pulse
% TDR_results.WC_ptdr_samples   worst case reflection samples of the reflection pulse
% TDR_results.ERL               reported effective return loss
% 
db = @(x) 20*log10(abs(x));
rms =@(x) norm(x)/sqrt(length(x));
if isfield(OP,'TDR_duration')
    TDR_duration=OP.TDR_duration; % approximate transit time multipler
else
    TDR_duration=5;
end
if ~isfield(OP,'DISPLAY_WINDOW')
    OP.DISPLAY_WINDOW=1; % approximate transit time multipler
end
f=sdd.Frequencies;
TDR_results.f=f;

% if license('test','RF_toolbox')
%     sdd.Parameters=s2s(sdd.Parameters,sdd.Impedance,ZT*2);
%     RL=squeeze(sdd.Parameters(param.RL_sel,param.RL_sel,:));
% else
for i = 1:length(sdd.Frequencies)
    rho= (2*ZT - sdd.Impedance) / (2*ZT + sdd.Impedance);
    interim =  sqrt(1-abs(rho)^2) * (1 - rho) / abs(1-rho);
    RL(i) = interim \ (sdd.Parameters(param.RL_sel,param.RL_sel,i) - rho) / ...
        (1 - rho *sdd.Parameters(param.RL_sel,param.RL_sel,i) ) * interim;
end
% end
RL=squeeze(RL);
f9=f/1e9;
tr=param.TR_TDR;
TDR_results.delay=500e-12 ;
% determine max time from thue
if sdd.NumPorts==1
    try
        maxtime=median(-unwrap(angle(squeeze(sdd.Parameters(1,1,:))))./sdd.Frequencies(:))*TDR_duration;
    catch
        maxtime=2e-9;
    end
    pix=1;
else
    [ fir4del, tu] =get_RAW_FIR(squeeze(sdd.Parameters(2,1,:)),f,OP,param);
    pix=find(fir4del==max(fir4del),1);
    maxtime=tu(pix)*TDR_duration+TDR_results.delay;
    if maxtime > tu(end); maxtime=tu(end);end
end
% add delay 500 ps for TDR and 3 times Gaussnan transtion time
% (makes gausian edge somewhat causal)

H_t = exp( -2*(pi*f9*(tr)/1.6832).^2 ).*exp(-(1j)*2*pi*f9*TDR_results.delay/1e-9).*exp(-1j*2*pi*f9*tr*3);
if ~isfield(OP,'cb_Guassian')
    Use_gaussian=1;
else
    Use_gaussian=OP.cb_Guassian;
end
if Use_gaussian
    if iscolumn(H_t), H_t=H_t.'; end
    RLf=RL(:).'.*H_t;
else % add extra 3x tr delay for causality
    RLf=RL.*exp(-(1j)*2*pi*f9*TDR_results.delay/1e-9);
end
if isfield(OP,'cb_BesselThompson')
    OP.BesselThompson=OP.cb_BesselThompson;
else
    OP.BesselThompson=0;
end
if isfield(param,'f_r')
    param.fb_BW_cutoff=param.f_r;
end
if OP.BesselThompson
    a = bessel( param.BTorder );
    acoef=fliplr( a );
    H_bt =a(1)./ polyval(acoef, (1i*f./(param.fb_BT_cutoff*param.fb)));
    %     acoef=[1    10    45   105   105];
    %     H_bt =acoef(5)./ polyval(acoef, (1i*f./(param.fb_BT_cutoff*param.fb)));
else
    H_bt=ones(1,length(f));
end

if isfield(OP,'Butterworth')
    if OP.Butterworth
        H_bw = 1./polyval([1 2.613126 3.414214 2.613126 1], 1i*f./(param.fb_BW_cutoff*param.fb));
    else
        H_bw=ones(1,length(f));
    end
else
    H_bw=ones(1,length(f));
end
if iscolumn(H_bt), H_bt=H_bt.';end
if iscolumn(H_bw), H_bw=H_bw.';end
if iscolumn(RLf), RLf=RLf.';end
TDR_results.Rx_filter=H_bt.*H_bw;
RLf=RLf.*TDR_results.Rx_filter;
TDR_results.tx_filter=H_t;


[IR, t, causality_correction_dB, truncation_dB] = ...
    s21_to_impulse_DC(RLf,  sdd.Frequencies(:), param.sample_dt,OP);


%
% param.tfx =4.2e-10; % need to put in xls file and comment this out
tfx=param.tfx;

 IR(abs(IR)<=OP.impulse_response_truncation_threshold)=0; % noise filter

t = t-TDR_results.delay;
tend=find(t>=maxtime,1);
if isempty(tend), tend=length(t); end
IR=IR(1:tend);
t=t(1:tend);
if isempty(tend), tend=length(t); end
tstart=find(t>=tr*1e-9,1); % account for Gaussian precursor
if isempty(tstart), tstart=1;end
if isempty(tend) || tstart >= tend  
    if isempty(tend) || tstart >= tend
        warndlg('TDR compuation not valid, try decreasing truncation tolerance, increasing samples, or adding a transmisson line','WrnTDR');
    end
    tend=length(t);
    tstart=1;
end
OP.cb_step=0; % step is a basically a cos^2 form -pi/4 to 0. not activated. 
ch=get_StepR(IR(tstart:tend),param,OP.cb_step,ZT);
TDR_results.tdr= ch.ZSR;
TDR_results.t = t(tstart:tend);
PTDR=get_PulseR(IR(tstart:tend),param,OP.cb_step,ZT);
if OP.TDR || OP.PTDR % determin average impededance with
    try
        tfstart=find(t>=3*tr*1e-9,1);
        x=squeeze(TDR_results.t(tfstart:end));TDR_results.x=TDR_results.x(:);
        y=squeeze(TDR_results.tdr(tfstart:end));TDR_results.y=TDR_results.y(:);
        w= exp(-(x-x(1))/OP.T_k ) ; % weighting function
        TDR_results.avgZport=mean(y.*w)/mean(w);
    catch
        TDR_results.avgZport=0;
        fit=zeros(1,1);
        p=[0 0 0 0  ];
    end
    TDR_results.RL=RL;
end
if OP.PTDR
    RL_equiv=-inf;
    L=param.levels;
    BinSize=OP.BinSize;
%     param.specBER=1e-5;
    if OP.DISPLAY_WINDOW
        hwaitbar=waitbar(0);
    else
        fprintf('Worst ERL searching');
    end
    % adjust PTDR for NDFE
    ntx=find(TDR_results.t >= tfx,1,'first');
    gatestartt=TDR_results.t(ntx);
    gatestartV=PTDR.pulse(ntx);
    ndfex=find(TDR_results.t > (param.ndfe+1)*param.ui+tfx,1,'first');
    tk=param.ui*1*(param.ndfe+1)+tfx;
    if  isempty(ndfex), ndfex=length(TDR_results.t); end
    PTDR.pulse_orig=PTDR.pulse;
    if(param.Grr_limit == 0 )
        fctrx(1:length(PTDR.pulse_orig))=(1+param.rho_x)*param.rho_x;
    else
        fctrx(1:length(PTDR.pulse_orig))=1;
    end
    if param.ndfe>0 && param.beta_x~=0;
        for ii=ntx:ndfex
            % adjust for near end loss
            fctrx0(ii)= 10.^(param.beta_x*(TDR_results.t(ii)-tk)/20);
            x=(TDR_results.t(ii)-tfx)/param.ui;
            % determine how much of the return loss to use base on expected
            % missing reflections
            fctrx1(ii)= (1+param.rho_x)*param.rho_x*exp(-(x-1*param.ndfe-1).^2/(1+param.ndfe)^2);
            fctrx(ii)=fctrx0(ii).*fctrx1(ii);
        end
        fctrx(1:ntx)=0;
        if isrow(fctrx), fctrx=fctrx(:);end
        PTDR.pulse=PTDR.pulse.*fctrx;
        if 0
            figure(10101)
            plot(TDR_results.t(1:ndfex)/param.ui,fctrx(1:ndfex),'disp','fctrx');
            hold on
            plot(TDR_results.t(1:ndfex)/param.ui,fctrx0(1:ndfex),'disp','fctrx0');
            plot(TDR_results.t(1:ndfex)/param.ui,fctrx1(1:ndfex),'disp','fctrx1');
        end
    end
    FAST_NOISE_CONV=0;
    ERLRMS=rms(PTDR.pulse);
    for ki=1:param.samples_per_ui
        progress = ki/param.samples_per_ui;
        if OP.DISPLAY_WINDOW
            waitbar(progress, hwaitbar, 'ERL computing'); figure(hwaitbar); drawnow;
        else
            if ~mod(progress*100,1), fprintf('%i%% ', progress*100 );end
        end
        tps=PTDR.pulse(ki:param.samples_per_ui:end);
        if OP.RL_norm_test
            rl_fom=(norm(tps));
        else
            testpdf=get_pdf_from_sampled_signal( tps,L, BinSize*10, FAST_NOISE_CONV );
            cdf_test=cumsum(testpdf.y);
            rl_test=(-testpdf.x(find(cdf_test>=param.specBER,1,'first')));
            rl_fom=rl_test;
        end
        if rl_fom > RL_equiv;
            RL_equiv=rl_fom;
            best_ki=ki;
        end
        if ~OP.RL_norm_test
            best_erl=rl_test;
            best_pdf=testpdf;
            best_cdf=cdf_test;
        end
        
    end
    if OP.RL_norm_test
        tps=PTDR.pulse(best_ki:param.samples_per_ui:end);
        testpdf=get_pdf_from_sampled_signal( tps,L, BinSize*10, FAST_NOISE_CONV );
        cdf_test=cumsum(testpdf.y);
        best_erl=(-testpdf.x(find(cdf_test>=param.specBER,1,'first')));
    end
    
    
    try
        close(hwaitbar)
    catch
    end
    if ~exist('best_ki','var'),best_ki=1;end
    TDR_results.ptdr_RL=PTDR.pulse;  % reflection waveform from the pulse
    TDR_results.WC_ptdr_samples_t=TDR_results.t(best_ki:param.samples_per_ui:end);
    TDR_results.WC_ptdr_samples=PTDR.pulse(best_ki:param.samples_per_ui:end);
    TDR_results.ERL=-db(best_erl);
    TDR_results.ERLRMS=-db(ERLRMS);
    
end


% end get TDR
%%



function result=get_StepR(ir,param,cb_step,ZT)
%ir = impulse response
%t_base=time array with equal time steps
%samp_UI = number of samples per UI for ir
%   result.SBR
% t for debug
t=(1/param.fb)/param.samples_per_ui*(0:length(ir)-1);

if cb_step
    Ag=1;
    dt=1/param.fb/param.samples_per_ui;
    edge_time=param.TR_TDR*1e-9;
    fedge=1/edge_time;
    tedge=0:dt:edge_time*2;
    %
    edge=Ag*(2*cos(2*pi*(tedge)*fedge/16-pi/4).^2-1);
    drive_pulse=[edge ones(1,param.samples_per_ui)];
    %pulse=filter(UI_ones,1,ir);
    
    pulse=filter(drive_pulse,1,ir);
else
    pulse=filter( ones(1,length(ir)),1,ir);
end
TDR_response=(1+pulse)./(1-pulse)*ZT*2;
result.ZSR=TDR_response;
result.pulse=pulse;


function result=get_PulseR(ir,param,cb_step,ZT)
%ir = impulse response
%t_base=time array with equal time steps
%samp_UI = number of samples per UI for ir

% t for debug
t=(1/param.fb)/param.samples_per_ui*(0:length(ir)-1);

if cb_step
    Ag=1;
    dt=1/param.fb/param.samples_per_ui;
    edge_time=param.TR_TDR*1e-9;
    fedge=1/edge_time;
    tedge=0:dt:edge_time*2;
    %
    edge=Ag*(2*cos(2*pi*(tedge)*fedge/16-pi/4).^2-1);
    drive_pulse=[edge ones(1,param.samples_per_ui)];
    %pulse=filter(UI_ones,1,ir);
    % t for debug
    t=(1/param.fb)/param.samples_per_ui*(0:length(ir)-1);
    
    pulse=filter(drive_pulse,1,ir);
else
    pulse=filter( ones(1,param.samples_per_ui),1,ir);
end
PDR_response=(1+pulse)./(1-pulse).*ZT*2;
result.PDR=PDR_response;
result.pulse=pulse;



function [ FIR t] =get_RAW_FIR(H,f,OP,param)
H_r = 1./polyval([1 2.613126 3.414214 2.613126 1], 1i*f./(0.75*param.fb));
if ~iscolumn(H), H=H.';end
if ~iscolumn(H_r), H_r=H_r.';end
H=H(:).*H_r;
[FIR, t, ~,~] = s21_to_impulse_DC(H ,f, param.sample_dt, OP) ;
% SBR=filter(ones(1, param.samples_per_ui), 1, FIR);

function [ a] = bessel( n )
% bessel polynomial
for ii= 0:n
    a(ii+1) = factorial(2*n-ii) /  (2^(n-ii)*factorial(ii)*factorial(n-ii));
end


function [ cmd_str ] = save_cmd_line( config_file, chdata, num_fext,num_next )
% save commmend string 
% for saving from interactive queries


cmd_str=[ mfilename '('  '''' config_file '''' ',' num2str(num_fext) ', ' num2str(num_next) ',' '''' chdata(1).filename ''''];
for i=1:num_next+num_fext
    cmd_str= [cmd_str  ',' '''' chdata(i+1).filename ''''];
end
cmd_str= [ cmd_str ')'];


%%%%% require the RF tool box
%%
function SLD=SL(S,f,R)
% source load impact return loss add to S21
% pass in S21
zref=100;
if R==0
    warndlg('Termination should not be set to zero');
        SLD=S;
    return
end
  
if R > zref
    spr =R_series2(zref,f,(R-zref)); % make series sparameter
    SLD=sparameters(cascadesparams(S.Parameters,spr.Parameters),f,S.Impedance); % Casdeade
elseif R < zref
    spr =r_parrelell2(zref,f,-R*zref/(R-zref)); % make parrellel sparameter
      SLD=sparameters(cascadesparams(S.Parameters,spr.Parameters),f,S.Impedance); 
else
    SLD=S;
end
%%

function S =R_series2(zref,f,R)
r=ones(1,length(f))*R;
S_Parameters(1,1,:) =  r./(r + 2*zref);
S_Parameters(2,2,:) =  r./(r + 2*zref);
S_Parameters(2,1,:) = (2*zref)./(r + 2*zref);
S_Parameters(1,2,:) =  (2*zref)./(r + 2*zref);
S=sparameters(S_Parameters,f,zref);

function S =r_parrelell2(zref,f,rpad)
S_Parameters(1,1,:) =  -zref/(rpad*(zref/rpad + 2)).*ones(1,length(f));
S_Parameters(2,2,:) =  -zref/(rpad*(zref/rpad + 2)).*ones(1,length(f));
S_Parameters(2,1,:) =   2/(zref/rpad + 2).*ones(1,length(f));
S_Parameters(1,2,:) =  2/(zref/rpad + 2).*ones(1,length(f));
S=sparameters(S_Parameters,f,zref);





function S =s_for_c2(zref,f,cpad)
% S is 2 port s parameters out
S_Parameters(1,1,:) =  -1i*2*pi.*f*cpad*zref./(2+1i*2*pi.*f*cpad*zref);
S_Parameters(2,2,:) =  -1i*2*pi.*f*cpad*zref./(2+1i*2*pi.*f*cpad*zref);
S_Parameters(2,1,:) = 2./(2+1i*2*pi.*f*cpad*zref);
S_Parameters(1,2,:) = 2./(2+1i*2*pi.*f*cpad*zref);
S=sparameters(S_Parameters,f,zref);

function S =s_for_c4(zref,f,cpad)

S2 = s_for_c2(zref,f,cpad);
S4P=s2_to_s4(S2.Parameters);
S=sparameters(S4P,f,zref);
S.Parameters=snp2smp(S.Parameters,zref,[ 1 3 2 4]);




%%
function [chdata, param] = process_sxp(param, OP, chdata, SDDch);    
num_files=length(chdata);
if ~OP.DISPLAY_WINDOW, fprintf('reading file '); end
for i=1:num_files
    if param.package_testcase_i==1  && i==1
         if OP.TDR && i==1
                S.Parameters(1,1,:)=SDDch(:,1,1);
                if ~param.FLAG.S2P;
                    S.Parameters(1,2,:)=SDDch(:,1,2);
                    S.Parameters(2,1,:)=SDDch(:,2,1);
                    S.Parameters(2,2,:)=SDDch(:,2,2);
                    S.Impedance=100;
                    S.NumPorts=4;
                else
                    S.Impedance=100;
                    S.NumPorts=1;
                end
                S.Frequencies=chdata(i).faxis;
                if OP.ADJ_4_LOAD
                    R_diepad = param.R_diepad;
                    C_diepad = param.C_diepad;
                    C_pkg_board = param.C_pkg_board;
                    pkg_Z_Tx=param.pkg_Z_c(1);
                    pkg_Z_Rx=param.pkg_Z_c(2);
                    % RX package length is assumed to be the same for all
                    % channel types.  swap Cp and Cd for Rx. ADJ_4_LOAD only
                    % for  Rx pkg
                    [ s11in, s12in, s21in, s22in ]= make_pkg(S.Frequencies, param.Pkg_len_RX,  C_pkg_board(2) ,C_diepad(2), pkg_Z_Rx, param);
                    S.Parameters(1,1,:)=s11in;
                    S.Parameters(1,2,:)=s12in;
                    S.Parameters(2,1,:)=s21in;
                    S.Parameters(2,2,:)=s22in;
                    S=sparameters(S.Parameters,S.Frequencies,100);
                    S=SL(S,S.Frequencies,R_diepad(2)*2);
                end
                
                % need to combine S wiht is page and channel
                if param.FLAG.S2P
                    port_sel=1;
                else
                    port_sel=[1 2];
                end
                OP.impulse_response_truncation_threshold=1e-5;
                for ipsl=1:length(port_sel) % do for both port if s4p
                    for izt=1:length(param.Z_t) % do for all tdr impedances
                        param.RL_sel=port_sel(ipsl);
                        OP.interp_sparam_phase='interp_to_DC'; % better for return loss
                        OP.interp_sparam_mag='trend_to_DC';
                        TDR_results(izt,ipsl) = get_TDR(S, OP, param,param.Z_t(izt));
                        if ipsl ==1
                            if  param.Z_t(izt) == 50
                                chdata(i).TDR11.ZSR=[TDR_results(izt,1).tdr];
                                chdata(i).TDR11.t=TDR_results(izt,1).t;
                            end
                            chdata(i).PDTR11(izt).ptdr=TDR_results(izt,1).ptdr_RL;
                        else
                            if  param.Z_t(izt) == 50
                                chdata(i).TDR22.ZSR=[TDR_results(izt,2).tdr];
                                chdata(i).TDR22.t=TDR_results(izt,2).t;
                            end
                            chdata(i).PDTR22(izt).ptdr=TDR_results(izt,2).ptdr_RL;
                                    end
                    end
                end
                chdata(i).TDR11.avgZport=[TDR_results(:,1).avgZport];
                if ~param.FLAG.S2P 
                    chdata(i).TDR22.avgZport=[TDR_results(:,2).avgZport];
                end
            if OP.PTDR && i==1
                chdata(i).TDR11.ERL=[TDR_results(:,1).ERL];
                chdata(i).TDR11.ERLRMS=[TDR_results(:,1).ERLRMS];
                if ~param.FLAG.S2P
                    chdata(i).TDR22.ERL=[TDR_results(:,2).ERL];
                    chdata(i).TDR22.ERLRMS=[TDR_results(:,2).ERLRMS];
                else
                    chdata(i).TDR22.ERL=[];
                    chdata(i).TDR22.ERLRMS=[];
                end
            else
                chdata(i).TDR11.ERL=[];
                chdata(i).TDR22.ERL=[];
                chdata(i).TDR11.ERLRMS=[];
                chdata(i).TDR22.ERLRMS=[];
            end
         end
 
    end
    if OP.DISPLAY_WINDOW && OP.DEBUG && OP.TDR
        h=figure(180);
        if param.package_testcase_i==1 && i == 1
            if i==1,
                htabgroup = uitabgroup(h);
                htab1 = uitab(htabgroup, 'Title', 'TDR TX');
                htab3 = uitab(htabgroup, 'Title', 'PTDR TX');
                hax1 = axes('Parent', htab1);
                hax3 = axes('Parent', htab3);
                if ~param.FLAG.S2P
                    htab2 = uitab(htabgroup, 'Title', 'TDR RX');
                    htab4 = uitab(htabgroup, 'Title', 'PTDR RX');
                    hax2 = axes('Parent', htab2);                   
                    hax4 = axes('Parent', htab4);                    
                end
            end
            set(h,'CurrentAxes',hax1)
            hold on
            plot(chdata(i).TDR11.t(:),chdata(i).TDR11.ZSR,'disp',[ chdata(i).base ' Tx port']);
            hold off
            legend (hax1, 'off');grid on;zoom xon;
            set(legend (hax1, 'show'), 'interp', 'none');

            if ~param.FLAG.S2P
                set(h,'CurrentAxes',hax2)
                hold on
                plot(chdata(i).TDR22.t(:),chdata(i).TDR22.ZSR,'disp',[ chdata(i).base ' Rx port']);
                hold off
                legend (hax2, 'off');grid on;zoom xon;
                set(legend (hax2, 'show'), 'interp', 'none');
            end
            
            set(h,'CurrentAxes',hax3)
            hold on
            for izt=1:length(param.Z_t)
                msg=[ chdata(i).base ' Tx port Zt='  num2str(param.Z_t(izt),3) '  ERL=', num2str(chdata(i).TDR11.ERL(izt), 3) 'db'];
                plot(chdata(i).TDR11.t,chdata(i).PDTR11(izt).ptdr,'disp',msg);
                msg=['PTDR Zt=' num2str(param.Z_t(izt),3) ' worst sampled noise cursors ' 'Tx port Zt='  num2str(param.Z_t(izt),3) ];
                stem(TDR_results(izt,1).WC_ptdr_samples_t,TDR_results(izt,1).WC_ptdr_samples,'disp',msg);
            end
            hold off
            legend (hax3, 'off');grid on;zoom xon;
            set(legend (hax3, 'show'), 'interp', 'none');
            if ~param.FLAG.S2P
                set(h,'CurrentAxes',hax4)
                hold on
                for izt=1:length(param.Z_t)
                    msg=[ chdata(i).base ' Tx port Zt='  num2str(param.Z_t(izt),3) '  ERL=', num2str(chdata(i).TDR22.ERL(izt), 3) 'db'];
                    plot(chdata(i).TDR22.t,chdata(i).PDTR22(izt).ptdr,'disp',msg);
                    msg=['PTDR Zt=' num2str(param.Z_t(izt),3) ' worst sampled noise cursors ' 'Rx port Zt='  num2str(param.Z_t(izt),3) ];
                    stem(TDR_results(izt,1).WC_ptdr_samples_t,TDR_results(izt,1).WC_ptdr_samples,'disp',msg);
                end
                hold off
                legend (hax4, 'off');grid on;zoom xon;
                set(legend (hax4, 'show'), 'interp', 'none');
            end
        end
    end
    if param.FLAG.S2P, return; end
end
