%% ===========================================================================
% Copyright (C):
%       2001-2002 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                         All rights reserved.
% Description   : ETSI ADSL performance simulations
%
% Project       : B1, FTW
% Author(s)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%
% Reference:
%    ETSI STC TM6 Draft RTS/TM-06006 (ADSL), rev 10a, 2001-05
%    ETSI STC TM6 Draft RTS/TM-06025 (ADSL), rev 1a, 2001-09
%    ETSI TM6 ADSL-EED4, which is available as part of TD34 at the ESTI TM6
%    meeting in Stockholm (013t34).
%    WD26R1 at the ESTI TM meeting in Stockholm (013w26r1).
%    WD20   at the ESTI TM meeting in Sophia Antipolis (014w20).
%    TD13   at the ESTI TM meeting in Sophia Antipolis (014t13).
%    ETSI TM6 021t31r1, Torino, Italy, 4th-8th Feb., 2002
%
% CVS:       $Id: ADSLperf.m,v 1.5 2002/03/13 08:19:50 tono Exp $
%% ===========================================================================

%% ===========================================================================
% Change History
%      2001-07-03 (ToNo) Created
%      2001-07-06 (ToNo) First time loop over all experiments was working
%      2001-07-19 (ToNo) Fixed the Zterm to 100 Ohm
%      2001-08-14 (ToNo) Fixed how to have both noise and same-pair-ISDN 
%      2001-08-17 (ToNo) Implemented a more robust reach search algorithm
%      2001-08-19 (ToNo) Fixed the minvalue for NT background noise
%      2001-09-05 (ToNo) Added loop dependent implementation loss
%      2001-09-12 (ToNo) Made a switch for dynamic BG noise
%      2001-09-12 (ToNo) Add support for "Maximum aggregate transmit power"
%      2001-12-14 (ToNo) Add nicer save and logging support
%      2001-12-17 (ToNo) Add support for noise injection normalization
%      2001-12-18 (ToNo) Separated background and modem noise
%      2002-01-11 (ToNo) Changed the way same pair ISDN is added
%      2002-01-21 (ToNo) Fixed the Zterm to 135 Ohm
%      2002-02-01 (ToNo) Added support for power cut-back
%      2002-02-06 (ToNo) Added a rough model of ISI/ICI
%% ===========================================================================
global ex;
global result;

%% ===========================================================================
% Logging, print, saving and plot options
format compact;
plotit   = 0;          % Set to one if you want to see some plots
maxplotf = 1.2e6;      % Max frequency to plot

dosingle = 0;          % Set to 1 if we want a single/smaller set of tests

% Get current date
dvec     = datevec(now);
isodate  = sprintf('%4d%02d%02d',dvec(1),dvec(2),dvec(3));

% Start logging
savefile = sprintf('matlab%s_perf',isodate);
logfile  = sprintf('diary%s_perf',isodate);
diary(logfile);

% Print greetings
fprintf('\nFTW ETSI ADSL performance simulations (%s)\n\n',...
        sprintf('%4d-%02d-%02d',dvec(1),dvec(2),dvec(3)));

%% ===========================================================================
% Get default structures for ETSI simulations
nomodems   = 0;                                % No of self disturbers (ADSL)
ex.param   = setupParam;                       % Get basic parameters
[ex.tfplist, ex.param.HAMBandName] = ...
	etsi_tfplanHAM([]);                    % Need a HAM band definition
ex.tfplist = fsan_tfplansMISC(ex.tfplist);     % Get plans for alien noise
ex.tfplist = etsi_noisesADSL(ex.tfplist);      % Get the ETSI ADSL noise models
ex.lclist  = setupLClist;                      % Line code defs (theoretical)
ex.lclist  = fsan_lcdefs(ex.lclist);           % Line code definitions (ADSL)
ex.clist   = etsi_cables([]);                  % Get ETSI (and ANSI) cables;
ex.tfplist = etsi_tfplansADSL(ex.tfplist);     % Get ETSI ADSL PSD tf plans
gui.ttlist = etsi_loopsADSL([]);               % Get ETSI ADSL loops

%% ========================================================================
% Connect tf-plans with names used in loop definitions

ISDNnoise = 'ETSI-ADSL-SamePairISDN4B3T';
isdntfplan = getList(ex.tfplist,ISDNnoise);

% Set up structure for same pair ISDN - will change during experiment
xDSL.name='SamePairISDN';  
xDSL.used='ETSI-ADSL-SamePairISDN4B3T';
ex.param.xDSLlist=insertList(ex.param.xDSLlist,xDSL);

noiseno=1;
noise_tfp(noiseno)=getList(ex.tfplist,xDSL.used);
noise_psd(noiseno).up=eval(noise_tfp(noiseno).PSD.upstream).*(10.^((ex.param.XTlevel.NEXT+15*log10(ex.param.frequency.f./1e6))/10));
noise_psd(noiseno).down=eval(noise_tfp(noiseno).PSD.downstream).*(10.^((ex.param.XTlevel.NEXT+15*log10(ex.param.frequency.f./1e6))/10));
noiseno=noiseno+1;

% Set up structure for ADSL - will change during experiment
xDSL.name='ADSL';  
xDSL.used='ADSL-ECoverPOTS';
tfp='ADSL-ECoverPOTS';
ex.param.xDSLlist=insertList(ex.param.xDSLlist,xDSL);

noise_tfp(noiseno)=getList(ex.tfplist,xDSL.used);
noise_psd(noiseno).up=eval(noise_tfp(noiseno).PSD.upstream);
noise_psd(noiseno).down=eval(noise_tfp(noiseno).PSD.downstream);
noiseno=noiseno+1;

% Set up structure for Noise - will change during experiment
xDSL.name='Noise';  
xDSL.used='ETSI-ADSL-ECoverPOTS Noise C'; %'TBD';
ex.param.xDSLlist=insertList(ex.param.xDSLlist,xDSL);

noise_tfp(noiseno)=getList(ex.tfplist,xDSL.used);
noise_psd(noiseno).up=eval(noise_tfp(noiseno).PSD.upstream);
noise_psd(noiseno).down=eval(noise_tfp(noiseno).PSD.downstream);
noiseno=noiseno+1;

if plotit,
figure;
f=ex.param.frequency.f;
semilogx(f,10*log10(noise_psd(2).up),...
     f,10*log10(noise_psd(2).down),...
     f,10*log10(noise_psd(3).up),...
     f,10*log10(noise_psd(3).down)...
     )
legend(sprintf('%s up',tfp), sprintf('%s down',tfp),...
       sprintf('%s up',xDSL.used), sprintf('%s down',xDSL.used));
title('PSD masks');
axis([1e3,maxplotf,-140,-20]) 
grid on;
drawnow;
end;

fprintf('Number of self disturbers = %d\n',nomodems);

%% ========================================================================
% Define some line code parameters
usdlcdef = 'ADSL-DMT';

framingoh                   = 0.032;  % Framing overhead of 32kbit/s
rsoh                        = 1.13;   % Read-Solomon overhead of 16%
tmp_lc=getList(ex.lclist,usdlcdef);
tmp_lc.param.xtalk_margin   = 6;      % We use a cross-talk margin of 6dB
tmp_lc.param.signal_margin  = 0;      % (and therefore no signal margin)
tmp_lc.param.SNRloss        = 1.5;    % SNR loss of 1.5dB - Changed later
tmp_lc.param.codingGain     = 4.25;   % Read-Solomon raw coding gain of 4.25dB
ex.lclist   = setList(ex.lclist,usdlcdef,tmp_lc); 
    
eval([tmp_lc.lcPrint,'(tmp_lc)']);    % Print line code setup
fprintf('Uses a %dkb/s framing overhead and a RS overhead of %g%%\n',framingoh*1e3,(rsoh-1)*100)

%% ========================================================================
% Decide what frequency axis to use
df		= 4.3125e3;           % Carrier spacing in Hz
if ex.param.frequency.fastrecalc,
    tmp=[];
else
    tmp=[df:df:256*df];
end;
ex.param.frequency.f = union(getMinFreqaxis(ex.param.frequency.max),tmp);
f=ex.param.frequency.f; 

%% ========================================================================
% Set up general parameters

ex.param.frequency.fastcalc  = 1;      % Use the fast frequency plan
ex.param.modemlist           = str2mat('ADSL'); % Modem(s) to be investigated
adslrespos                   = 1;      % ADSL result position in result struct
ex.param.frequency.max       = 1.2e6;  % Max frequency we want to use
ex.param.XTlevel.NEXTx4      = 1;      % Flag to use forth order next modeling
ex.param.bgNoise.LT          = '-140'; % Backgound noise at LT (dBm/Hz)
ex.param.bgNoise.NT          = '-140'; % Backgound noise at NT (dBm/Hz)
ex.param.modemNoise.LT       = '-120'; % Modem noise at LT (dBm/Hz)
ex.param.modemNoise.NT       = '-120'; % Modem noise at NT (dBm/Hz)
ex.param.Zterm               = 135;    % ADSL is using 100 Ohm but not for s21!

DoPowerCutBack               = 1;      % Adjust for power cut-back
DoAdjustForISCI              = 1;      % Adjust for ISI/ICI
DoAdjustForNID               = 1;      % Adjust for noise injection dev
DoAdjustForMN                = 1;      % Adjust for modem noise
useDynamicMN                 = 1;      % Use the dynamic modem noise
                                       % (loop length and frequency dependent)

% Loop definition contain an offset and a position of the variable segment
offsetlen=[0 0 1500 2200 1750 1750 4200 1100];
ypos=[ 1 1 1 2 1 1 2 1]+1; % position of the variable tt.topology segment

if DoAdjustForNID,
    fprintf('With adjustment for noise injection normalization\n')
end;

if DoAdjustForMN,
    fprintf('With a modem noise floor\n')
end;

if DoPowerCutBack,
    fprintf('With a power cut back\n')
end;

if DoAdjustForISCI,
    fprintf('With addition of ISI/ICI noise\n')
end;

%% ===========================================================================
% Set up experiment parameters and loops

% Start timer
tic;
noiters=0;
totalcnt=0;
nonoreach=0;

if dosingle,
    ecs=0;
    pots=0;
    noises=['A'];
    direct=[0];
    loops=8; %loops=1:8;
    singlebitrate=128/1e3;%6144/1e3;
else % All experiments, we have 2*2*4*8*(6+10)=2048 reaches to calculate
    ecs=[1 0];                         % 1/0 with or without echo cancellation
    pots=[0 1];                        % 1/0 over or not over POTS
    noises=['A','B','C','D'];          % Noise models to use (A-D)
    direct=[1 0];                      % 1/0 test DS or US direction
    loops=1:8;                         % Loops to test (1-8)
    singlebitrate=0;                   % 0 if multiple (all) bitrates is used
end

%% ===========================================================================
% Begin the nested loops

for withEC = ecs

if withEC,
    withNewAssumptions = 0; 
    DoAdjustForISCI_FDDoverISDN = 0;
else
    % Use new assumption decided at ETSI Stockholm 2001 meeting (for FDD)
    withNewAssumptions = 1; 
end;

for overPOTS = pots    
for noise=noises
for checkDS=direct

fprintf('\n=========================\n');

if checkDS,
    dir='DS';
    bitrates=[128 512 768 1024 1544 2048 3072 4096 5120 6144]/1e3; % DS bitrates
else
    dir='US';
    bitrates=[64 128 256 384 512 640]/1e3; % US bitrates
end;
if singlebitrate~=0
    bitrates=singlebitrate;
end;

% SNR loss depending on loop
if withNewAssumptions,
    % FDD
    SNRlossListDS = [2 2 2.5 2.5 2.5 2.5 2.5 3];
    if overPOTS,
        SNRlossListUS = 2*ones(1,8);
        isci.maxSNR = 42;
        isci.lowdeltaSNR = 6;
        isci.highdeltaSNR = 9;
    else
        SNRlossListUS = [2 2 2 2.5 2 2 2 2.5];
        isci.maxSNR = 41;
        isci.lowdeltaSNR = 8;
        isci.highdeltaSNR = 8;
    end;
else
    % EC
    SNRlossListDS = 1.5*ones(1,8);
    SNRlossListUS = 1.5*ones(1,8);
end;

tabdata = zeros(length(loops),length(bitrates));

for targetrate=bitrates
for loopno=loops

%% ===========================================================================
% Use the right tfplan and set up allocation of tones

if withEC
    splitter='';
    if overPOTS
        tfp='ADSL-ECoverPOTS';
        tonesup=[7:31];
        tonesdown=[7:63,65:255];
        PxDown=20.4;
        PxUp=12.5;
    else
        tfp='ADSL-ECoverISDN';
        tonesup=[33:63];
        tonesdown=[33:95,97:255];
        PxDown=19.9;
        PxUp=13.3;
    end
else
   if withNewAssumptions,
       splitter='_SP';
   else
       splitter='';
   end;
   if overPOTS
        tfp='ADSL-FDDoverPOTS';
        tonesup=[7:30];
        tonesdown=[38:63,65:255];
        PxDown=19.9;
        PxUp=12.5;
        if DoAdjustForISCI,
            DoAdjustForISCI_FDDoverISDN=1;
        end;
    else
        tfp='ADSL-FDDoverISDN';
        tonesup=[33:56];
        tonesdown=[64:95,97:255];
        PxDown=19.3;
        PxUp=13.3;
        if DoAdjustForISCI,
            DoAdjustForISCI_FDDoverISDN=1;
        end;
    end
end    

% Change the lc according to the defined tones; SNR loss dep on loop
tmp_lc=getList(ex.lclist,usdlcdef);
tmp_lc.param.dmt.tonesus    = tonesup;
tmp_lc.param.dmt.tonesds    = tonesdown;

tmp_lc.param.PxDown         = PxDown;
tmp_lc.param.PxUp           = PxUp;

if checkDS, % Set up correct SNR loss value
  tmp_lc.param.SNRloss      = SNRlossListDS(loopno);
else
  tmp_lc.param.SNRloss      = SNRlossListUS(loopno);
end;

ex.lclist   = setList(ex.lclist,usdlcdef,tmp_lc); 

% Set up correct frequency for the IL values
if checkDS,
  f_T = 300e3;
else
  if overPOTS
      f_T = 75e3;
  else
      f_T = 150e3;
  end;
end;

% Set up same pair noise (NONE or 4B3T ISDN)
if overPOTS
    samepairnoise='ETSI-ADSL-NoSamePairISDN';
else
    samepairnoise='ETSI-ADSL-SamePairISDN4B3T'; 
end

% Redefine the changed structures
ex.param.FSANNoiseModel = sprintf('ETSI-%s Noise %s',tfp,noise);
scenario = sprintf('ETSI-ADSL-Loop%d',loopno);
ex.tt    = getList(gui.ttlist,scenario);     % Get our example scenario
xDSL.name='SamePairISDN';  
xDSL.used=samepairnoise;
ex.param.xDSLlist=setList(ex.param.xDSLlist,xDSL.name,xDSL);
xDSL.name='ADSL';  
xDSL.used=[tfp splitter];
ex.param.xDSLlist=setList(ex.param.xDSLlist,xDSL.name,xDSL);
xDSL.name='Noise';  
xDSL.used=sprintf('ETSI-%s Noise %s',tfp,noise);
ex.param.xDSLlist=setList(ex.param.xDSLlist,xDSL.name,xDSL);

noisepos= length(ex.param.xDSLlist);

%% ========================================================================
% Print experiment setup
fprintf('\n');
fprintf('Test %s using Noise "%s" with %s target of %gkbit/s\n',...
        [tfp splitter],xDSL.used,dir,targetrate*1e3);

% Show traffic and topology structure as plots (if wanted)
if plotit,
figure; plotTTstructure(ex.tt);
tfplan=getList(ex.tfplist,[tfp splitter]);
fax.min=1e3;fax.max=1.2e6;
figure; plotTFplanDMT(tfplan,'Linear',fax,tonesup,tonesdown);
drawnow; % Show it now
end;

%% ========================================================================
% Prepare the search, by setting a lower and upper bound on the reach value

minlen=2;

%%%% First get the lower bound (just a minimum distance)
xlen=minlen;    
ex.tt.topology{ypos(loopno),1}=xlen;
% Modem noise at NT depend on loop length 
if useDynamicMN,
    if withNewAssumptions,
        ex.param.modemNoise.NT = ...
            sprintf('get_noise_floor(ex.param.frequency.f, %f, %d)',...
            xlen+offsetlen(loopno),~overPOTS);
    else    
        ex.param.modemNoise.NT = ...
            sprintf('%g',max(-140,-20*(xlen+offsetlen(loopno))/4000-120));
    end;
end;

%%% Evaluate experiment %%%
result = evalExperiment;
% Adjust for ADSL performance tests
if ~overPOTS,      AddNoisePSD(isdntfplan); end  % Add same pair ISDN
if DoPowerCutBack, AdjustForPOB; end;            % Adjust for power cut back
if DoAdjustForNID, AdjustForNID; end;            % Adjust for noise inject.dev.
if DoAdjustForISCI_FDDoverISDN,AdjustForISCI; end; % Adjust for ISI/ICI
if DoAdjustForMN,  AdjustForMN; end;             % Floor for modem noise
% Get the bitrate
[bitrate_LT, bitrate_NT] = calcFSANresult(ex,result);
if checkDS,
    bitrate=raw2user(bitrate_NT(adslrespos));
else
    bitrate=raw2user(bitrate_LT(adslrespos));
end;
cnt=1;
low_xlen=xlen;
low_rate=bitrate;

%%% Now get the upper bound
if bitrate < targetrate
    % If we do not get the target rate at minimal reach then we are done
    quitloop=1;
    nonoreach=nonoreach+1;
    high_xlen=low_xlen-1;
    high_rate=low_rate-1;
else
    % Otherwise get a second distance as an upper bound
    quitloop=0;
    xlen=3000;
    contsearch=1;
    while contsearch
        % Modem noise at NT depend on loop length 
        if useDynamicMN,
            if withNewAssumptions,
                ex.param.modemNoise.NT = ...
                   sprintf('get_noise_floor(ex.param.frequency.f, %f, %d)',...
                   xlen+offsetlen(loopno),~overPOTS);
            else    
                ex.param.modemNoise.NT = ...
                   sprintf('%g',max(-140,-20*(xlen+offsetlen(loopno))/4000-120)); 
            end;    
        end;
        ex.tt.topology{ypos(loopno),1}=xlen; 
        result = evalExperiment;
        % Add same pair ISDN:
        if ~overPOTS,      AddNoisePSD(isdntfplan); end  
        % Adjust for power cut back
        if DoPowerCutBack, AdjustForPOB; end;
        % Adjust for noise normalization mismatch:
        if DoAdjustForNID, AdjustForNID; end;
        % Adjust for ISI/ICI
        if DoAdjustForISCI_FDDoverISDN,AdjustForISCI; end; 
        % Adjust for modem noise:
        if DoAdjustForMN,  AdjustForMN; end;
        [bitrate_LT, bitrate_NT] = calcFSANresult(ex,result);
        if checkDS,
            bitrate=raw2user(bitrate_NT(adslrespos));
        else
            bitrate=raw2user(bitrate_LT(adslrespos));
        end; 
        cnt=cnt+1;
        if bitrate<=0,
            % To far out -- no bitrate at all
            xlen = low_xlen + (xlen-low_xlen)/2;
        elseif bitrate<targetrate
            % A useful distance outside reach, i.e done in this phase
            contsearch=0;
            high_rate=bitrate;
            high_xlen=xlen;
        else
            % Still above target bitrate aim higher and move minimum distance
            lenrate=(low_xlen-xlen)/(low_rate-bitrate);
            low_rate=bitrate;
            low_xlen=xlen;
            xlen=low_xlen + lenrate*(targetrate-bitrate) + 100;
        end;
    end;
end;

%% ========================================================================
% Now search for the reach, between upper and lower bound, 
% downto a 1 m accuracy.

xlen=low_xlen + (high_xlen-low_xlen)/2;
while ~quitloop
    cnt=cnt+1;
    ex.tt.topology{ypos(loopno),1}=xlen;
    % Modem noise at NT depend on loop length 
    if useDynamicMN,
        if withNewAssumptions,
            ex.param.modemNoise.NT = ...
                sprintf('get_noise_floor(ex.param.frequency.f, %f, %d)',...
                xlen+offsetlen(loopno),~overPOTS);
        else    
            ex.param.modemNoise.NT = ...
                sprintf('%g',max(-140,-20*(xlen+offsetlen(loopno))/4000-120));
        end;
    end;
    result = evalExperiment;
    if ~overPOTS,      AddNoisePSD(isdntfplan); end  % Add same pair ISDN
    if DoPowerCutBack, AdjustForPOB; end;
    if DoAdjustForNID, AdjustForNID; end;
    if DoAdjustForISCI_FDDoverISDN,AdjustForISCI; end; % Adjust for ISI/ICI
    if DoAdjustForMN,  AdjustForMN; end;
    [bitrate_LT, bitrate_NT] = calcFSANresult(ex,result);
    if checkDS,
        bitrate=raw2user(bitrate_NT(adslrespos));
    else
        bitrate=raw2user(bitrate_LT(adslrespos));
    end;

    % Did we get a new lower or upper bound?
    if bitrate>targetrate,
        low_rate=bitrate;
        low_xlen=xlen;
    else
        high_rate=bitrate;
        high_xlen=xlen;
    end;
    
    % Do a "binary search" step (next length is in the middle of high and low)
    xlen=low_xlen + (high_xlen-low_xlen)/2;
    if high_xlen-xlen < 0.5,
        quitloop=1;
    end;
end %while  

low_xlen=floor(low_xlen); % Floor it down to the closest meter

%% ========================================================================
% Print result
if low_xlen>minlen,
    fprintf('Loop %d, Len=%.0fm, IL = %.1fdB (@%.0fkHz)\n',...
            loopno,low_xlen+offsetlen(loopno),...
            ADSLIL(loopno,f_T,low_xlen+offsetlen(loopno)),f_T/1e3);
    % Save the data in a structure for later processing
    tabdata(loopno,find(bitrates==targetrate))=low_xlen+offsetlen(loopno);
else
    fprintf('Loop %d, - \n',loopno);
    % Save the data in a structure for later processing
    tabdata(loopno,find(bitrates==targetrate))=0;
end;

% Save the data in a structure for later processing
tabdata(loopno,find(bitrates==targetrate))=low_xlen+offsetlen(loopno);

%% ========================================================================
% Loop around
noiters=noiters+1;
totalcnt=totalcnt+cnt;

end %loop
end %rate
if overPOTS,
    atype='POTS';
else
    atype='ISDN';
end;
if withEC,
    etype='EC';
else
    etype='FDD';
end;

% Save the table of data in a variable for later processing
eval(sprintf('ETSI_ADSL_%s_%s_%s_%s=tabdata;',atype,etype,dir,noise));

end %DS/US
end %noise
end %POTS
end %EC

%% ========================================================================
% Print elapsed time
usedtime=toc; 
fprintf('\n');
fprintf('Used time = %.2fs (%.2fs per reach)\n',usedtime,usedtime/noiters);
fprintf('Calculated %d experiments ',totalcnt);
fprintf('(in avg. %.1f experiments per nonzero reach)\n',...
        totalcnt/(noiters-nonoreach));

%% ========================================================================
% Save and stop logging
save(savefile);
diary off;