%% ===========================================================================
% Copyright (C) 1999 by Telia Research AB, Lulea, Sweden; All rights reserved. 
% 
% Description   : Test of Mixed services scenarios
% Project       : FSAN xDSL simulation tool
% Author(s)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%
% CVS:       $Id: MixedServices.m,v 1.4 2000/05/26 13:32:44 uvan Exp $
%% ===========================================================================

%% ===========================================================================
% Change History
%      1999-09-29 (ToNo) Created from ETSI reach table calcs
%      1999-11-17 (ToNo) Cleaning up so this script can be released as an example
%      2000-02-07 (DaB)  debugging
%      2000-04-05 (UvAn) Added frequency axis for the call to calcPSD
%      2000-05-10 (UvAn) Added tmp_tfplan.timeDivision.sync = 1 variable
%      2000-05-26 (UvAn) Added possibillities to fix the bitrate 
%% ===========================================================================

%% ===========================================================================
%  In the ETSI contribution on mixed services (994tXXa0) I used the following
%  parameters:
%  Table 1: mss=2; planselect=45;  % That is, TEx with 4-band plan
%  Table 2: mss=2; planselect=783; % That is, TEx with 7-band plan
%  Table 3: mss=1; planselect=45;  % That is, TCab with 4-band plan
%  Table 4: mss=1; planselect=783; % That is, TCab with 7-band plan
%% ===========================================================================
mss        = 2;	 % Mixed service scenario to do, 1=TCab, 2=TEx
planselect = 45; % 45 is infineon 4-band, 783 is my scalable 7-band plan;

global ex;
global simures;

% Set up default experiment structures (as suggested by ETSI)
ETSI_VDSL_defs;

% Setup some default line-code parameter values
lcname = 'VDSL-theo';
lc=getList(ex.lclist,lcname);
%lc.param.shannonGap     = 9.8+6-3.8;    % Shannon gap + Noise margin - Coding gain (dB) 
lc.param.codingGain=3.8;
lc.param.xtalk_margin   = 6;
lc.param.refSNR         =9.8;

lc.param.efficiencyLoss = 0.1;          % Efficiency loss (0-1)
lc.param.SNRloss	= 0;
ex.lclist=setList(ex.lclist,lc.name,lc);

% Only test VDSL
ex.param.modemlist=['VDSL'];

% Experiments to conduct
loops=[1];			        % Loops to test
switch mss
    case 1, % Cab scenario, TCab
        Gdists=[700 500 1400 1000 800];   	% Modems at these distances
        Grates=[Srates(:,3:4),Arates(:,2:4)]';  % Goal rates
        Gnames=[Snames(3:4,:); Anames(2:4,:)];  % Service names
        FSANnoise={'A'};                        % Noise used is A
        psdset=2;                               % PSD set used is PCabDM2
        pbolens = [800 900 1000];               % Typical lengths to test
        pbofreqs = [1.2e6 1.7e6 2.8e6];         % Typical frequencies to test
    case 2, % Exchange scenario, TEx
        Gdists=[1250 1100 700 1500 1200];       % Modems at these distances
        Grates=[Srates(:,1:3),Arates(:,2:3)]';  % Goal rates
        Gnames=[Snames(1:3,:); Anames(2:3,:)];  % Service names
        FSANnoise={'E'};                        % Noise used is E
        psdset=1;                               % PSD set used is PExD1M1
        pbolens = [1000 1100 1200];             % Typical lengths to test
        pbofreqs = [1.2e6 1.7e6 2.8e6];         % Typical frequencies to test
end;

[dists,distorder]=sort(Gdists); % Modems needs to be in length order
nomodemspernode=4;		% Number of VDSL modems per node

% Set up the different PBO methods to test (more after plan is set)
% Test without PBO
pbomethods    = {'None'};
pboparam.len  = {[0]};
pboparam.freq = {[0]};
pboparam.maxlen=max(dists);
% Add tests with reference freqency (constant PBO)
for i=1:length(pbofreqs),
    pbomethods = {pbomethods{:} 'RefFreq'};
    pboparam.len  = {pboparam.len{:}  0};
    pboparam.freq = {pboparam.freq{:} [pbofreqs(i)]};
end;
% Add tests with reference length
for i=1:length(pbolens),
    pbomethods = {pbomethods{:} 'RefLen'};
    pboparam.len  = {pboparam.len{:}  [pbolens(i)]};
    pboparam.freq = {pboparam.freq{:} 0};
end;
% Add tests with reference FEXT
for i=1:length(pbolens),
    pbomethods = {pbomethods{:} 'RefFEXT'};
    pboparam.len  = {pboparam.len{:}  [pbolens(i)]};
    pboparam.freq = {pboparam.freq{:} 0};
end;
% Add a test with reference noise
pbomethods    = {pbomethods{:} 'RefNoise'};
pboparam.len  = {pboparam.len{:}  0};
pboparam.freq = {pboparam.freq{:} 0};


% PSD mask to be contain within
mask.base = 'VDSL P';
mask.ext.cab = 'CabD';
mask.ext.ex = 'ExD1';

% Control plots and other outputs
doplots=1;
printpower=1;

% Frequency plan settings
freqplanbase='NBand'; 

stepfq=.2e6; 	
stopstepfq=.05e6;

switch planselect
   case 201,
        % 2 band all frequencies at HAM bands high, based on 402
        initialfq= [.138e6 3.65e6 7.1e6];
        movablefq=[ ]; % Do not try to move anything
        fqdir='du';
    case 202,
        % 2 band all frequencies at HAM bands high, based on 402
        initialfq= [7.1e6 10.1e6  20e6];
        movablefq=[ ]; % Do not try to move anything
        fqdir='du';
    case 2782,
        % 2 band based on the 7 band plan #782
        initialfq= [8.8e6 12.65e6  20e6]; 
        fqdir='du';
        movablefq=[length(initialfq)-1]; 
    case 3,
        % 3 band below 4.4
        initialfq= [30e3 .138e6 2.45e6 4.4e6];
        movablefq=[ 3]; % only change these frequencies among the fq frequencies
        fqdir='udu';
    case 34,
        % 4 band below 4.4 (including ADSL upsteam);
        initialfq= [30e3 .138e6 1.104e6 2.0e6 4.4e6];
        movablefq=length(initialfq)-1;
        fqdir='udud';
    case 301,
        % 3 band all frequencies at HAM bands high reduced to 10MHz, based on 402
        initialfq= [.138e6 3.65e6 7.1e6 10e6];
        movablefq=[ ]; % Do not try to move anything
        fqdir='dud';
    
    case 4,
        % 4 band Best over all
        initialfq= [.138e6 2.6e6 5.4e6 10.25e6  20e6];
        movablefq=[ 2 3 4]; % only change these frequencies among the fq frequencies
        fqdir='dudu';
    case 417, 
        % 4 band Use only up to 17.66 from the best over all
        initialfq= [.138e6 2.6e6 5.4e6 10e6 17.66e6];
        movablefq=[]; % Do not try to move anything
        fqdir='dudu';
    case 411,
        % 4 band with fixation for all but the middle one
        initialfq= [.138e6 1.9e6 4.2e6 10.1e6  20e6];
        movablefq=[3]; % fix all but the middle one
        fqdir='dudu';
    case 412,
        % 4 band with fixation for all but the middle one
        initialfq= [.138e6 3.65e6 6.9e6 10.1e6  20e6];
        movablefq=[3]; % fix all but the middle one
        fqdir='dudu';
    case 401,
        % 4 band Fix all frequencies at HAM bands low
        initialfq= [.138e6 1.9e6 3.65e6 10.1e6  20e6];
        movablefq=[ ]; % Do not try to move anything
        fqdir='dudu';
    case 402,
        % 4 band Fix all frequencies at HAM bands high
        initialfq= [.138e6 3.65e6 7.1e6 10.1e6  20e6];
        movablefq=[ ]; % Do not try to move anything
        fqdir='dudu';
    case 45,
        % 4 band based on NT119
        initialfq= [.138e6 3.75e6 8.325e6 12.5e6 20e6]; 
        fqdir='dudu';
        movablefq=[2:length(initialfq)-1];         
    case 455,
        % 4 band based on NT119 + ADSL upstream
        initialfq= [30e3 .138e6 3.75e6 8.325e6 12.5e6 20e6]; 
        fqdir='ududu';
        movablefq=[2:length(initialfq)-1];         
    case 451,
        % 2 band based on NT119
        initialfq= [.138e6 3.75e6 8.325e6 ]; 
        fqdir='du';
        movablefq=[2:length(initialfq)-1];         
    case 452,
        % 2 band based on NT119
        initialfq= [8.325e6 12.5e6 20e6]; 
        fqdir='du';
        movablefq=[2:length(initialfq)-1];         
    case 453,
        % 3 band based on NT119
        initialfq= [.138e6 3.75e6 8.325e6 10e6]; 
        fqdir='dud';
        movablefq=[2:length(initialfq)-1];         
    case 47821,
        % 4 band based on the 7 band plan #783
        initialfq= [.138e6 2.45e6 4.4e6 6.3e6 8.8e6]; 
        fqdir='dudu';
        movablefq=[length(initialfq)-1]; 
    case 47822,
        % 4 band based on the 7 band plan #783
        initialfq= [4.4e6 6.3e6 8.8e6 12.65e6  17.66e6]; 
        fqdir='dudu';
        movablefq=[length(initialfq)-1]; 
        
    case 48,
        % 4 band below 8.8 excluding ADSL upstream
        initialfq= [.138e6 2.1e6 4.15e6 6.7e6  8.8e6];
        movablefq=[ 2 3 4]; % only change these frequencies among the fq frequencies
        fqdir='dudu';
    case 49,
        % 4 band Suggested by Alcatel
        initialfq= [.138e6 3.1e6 5.5e6 9.5e6 20e6];
        movablefq=[ 2 3 4]; % only change these frequencies among the fq frequencies
        fqdir='dudu';

    case 51,
        % 5 band (4 band plus ADSL upstream)
        initialfq= [30e3 .138e6 2.6e6  5.4e6 10.0e6  20e6]; 
        fqdir='ududu';
        movablefq= 3:(length(initialfq)-1); % Do not move the first two frequencies
    case 511,
        % 5 band (4 band plus ADSL upstream)
        initialfq= [30e3 .138e6 2.55e6  5.4e6 10.4e6  20e6]; 
        fqdir='ududu';
        movablefq= 3:(length(initialfq)-1); % Do not move the first two frequencies
    case 52,
        % 5 band (excluding ADSL upstream)
        initialfq= [.138e6 2.8e6  4.6e6 5.8e6 11.9e6  20e6]; 
        fqdir='dudud';
        movablefq= 2:(length(initialfq)-1); % Do not move the first two frequencies
    case 581,
        % 5 band below 8.8MHz
        initialfq= [30e3 .138e6 2.1e6 4.1e6 6.7e6 8.8e6];
        movablefq=3:(length(initialfq)-1);
        fqdir='ududu';
    case 582,
        % 5 band below 8.8MHz based on the 3 band solution
        initialfq= [30e3 .138e6 2.45e6 4.4e6 6.3e6 8.8e6];
        movablefq=length(initialfq)-1; % only move the last grid point
        fqdir='ududu';

    case 6,
        % 6 band (5 band plus ADSL upstream)
        initialfq= [30e3 .138e6 2.0e6  3.8e6 5.5e6 10.1e6  20e6]; 
        fqdir='ududud';
        movablefq= 3:(length(initialfq)-1); % Do not move the first two band
    case 681,
        % 6 band (5 band plus ADSL upstream)
        initialfq= [30e3 .138e6  1.104e6 2.75e6 4.4e6 6.3e6 8.8e6]; 
        fqdir='ududud';
        movablefq= length(initialfq)-1; 

    case 71,
        % 7 band
        %initialfq= [30e3 .138e6 2.0e6 4.3e6 7.1e6 10.1e6 14.4e6 20e6]; % Telia start
        initialfq= [30e3 .138e6 2.5e6 4.9e6 7e6 10.15e6 14.45e6 20e6]; % Best so far
        fqdir='udududu';
        movablefq= 3:(length(initialfq)-1); % Do not move the first two band
    case 72,
        % 7 band
        %initialfq= [30e3 .138e6 1.104e6 1.9e6 3.8e6 5.5e6 10e6  20e6]; % Alcatel start
        initialfq= [30e3 .138e6 1.7e6 1.9e6 3e6 5.7e6 11e6 20e6]; % Result
        fqdir='udududu';
        movablefq= 3:(length(initialfq)-1); % Do not move the first two band
    case 78,
        % 7 band based on the 5 band below 8.8MHz
        initialfq= [30e3 .138e6 2.6e6 5.1e6 6.85e6 8.8e6 12.95e6  20e6]; 
        fqdir='udududu';
        movablefq=[length(initialfq)-1]; % Do not move anything but the last
    case 781,
        % 7 band based on the 5 band below 8.8MHz (best in show so far!)
        initialfq= [30e3 .138e6 2.6e6 5.1e6 6.85e6 8.8e6 12.95e6  20e6]; 
        fqdir='udududu';
        movablefq=[3:length(initialfq)-1]; 
    case 782,
        % 7 band based on the 5 band below 8.8MHz
        initialfq= [30e3 .138e6 2.45e6 4.4e6 6.3e6 8.8e6 12.65e6 20e6]; 
        fqdir='udududu';
        movablefq=[length(initialfq)-1]; 
    case 783,
        % 7 band based on the 5 band below 8.8MHz
        initialfq= [30e3 .138e6 2.45e6 4.4e6 6.3e6 8.8e6 12.65e6 17.66e6]; 
        fqdir='udududu';
        movablefq=[length(initialfq)-1]; 
    
    case 881,
        % 8 band (7 band plus ADSL upstream) based on 681
        initialfq= [30e3 .138e6  1.104e6 2.75e6 4.4e6 6.3e6 8.8e6 13e6 20e6]; 
        fqdir='udududud';
        movablefq= length(initialfq)-1; 
    case 91,
        initialfq = [30e3 .138e6 1.104e6 1.5e6 1.9e6 2.52e6 3.13e6 3.75e6 ... 
                4.62e6 5.5e6 7.15e6 8.325e6 9.5e6 ...
                10.12e6 12.46e6 14.77e6 17.66e6];
        fqdir = 'uduuddduudddduuu'; % Global at ETSI
        movablefq=[];
    otherwise error('Unknown planselect');
end;

% Print experiment setup
fprintf('\nSimulating with:\n');
fprintf('Plan #%d\n',planselect);
eval([lc.lcPrint,'(lc)']); % Print line code setup

% Set up a default PSD mask for VDSL (and the change only PSD def later)
% Default values
tmp_tfplan.name                  = 'VDSL ETSI';
tmp_tfplan.PSD.downstream        = 'N/A';
tmp_tfplan.PSD.upstream          = 'N/A';
tmp_tfplan.PSD.active.upstream   = [30e3 30e6];
tmp_tfplan.PSD.active.downstream = [30e3 30e6];
tmp_tfplan.PSD.PBO.method        = 'None';	
tmp_tfplan.PSD.PBO.param.len     = 1500;	
tmp_tfplan.PSD.PBO.param.freq    = 2e6;	
tmp_tfplan.PSD.PBO.param.maxlen  = pboparam.maxlen;
tmp_tfplan.PSD.PBO.param.plan    = planConvertion(initialfq,fqdir,0);
tmp_tfplan.PSD.HAM.active        = 0;
tmp_tfplan.timeDivision.up       = 1; % Time used in up resp. down lin
tmp_tfplan.timeDivision.down     = 1;
tmp_tfplan.timeDivision.sync     = 1;
tmp_tfplan.fixBitrate.name       = 'None';  % Fixed bitrate method
tmp_tfplan.fixBitrate.active     = 0;       % 0=inactive, 1=active
tmp_tfplan.fixBitrate.param      = 'N/A';   
tmp_tfplan.lcname                = lcname;

ex.tfplist=insertList(ex.tfplist,tmp_tfplan);

% Set up test with multiband reference length
noupbands=length(tmp_tfplan.PSD.PBO.param.plan.up.fc);
maxl=max(Gdists);
minl=min(Gdists);
deltal=(maxl-minl)/noupbands;
lens=(maxl-deltal/2):-deltal:(minl+deltal/2);
pbomethods    = {pbomethods{:}    'RefLenM'};
pboparam.len  = {pboparam.len{:}  lens};
pboparam.freq = {pboparam.freq{:} zeros(1,noupbands)};
% Add tests with multi reference length when the above "algorithm" isn't
% doing to well setting the lengths
if planselect==45,
pbomethods = {pbomethods{:} 'RefLenM'};
pboparam.len  = {pboparam.len{:}  [1100 800]};
pboparam.freq = {pboparam.freq{:} [0 0]};
end;

% Choose VDSL scheme to be evaluated
vdslDuplex = tmp_tfplan.name;
tmp_xDSL.name='VDSL';
tmp_xDSL.used=vdslDuplex;
ex.param.xDSLlist(1)=tmp_xDSL; 

fq=initialfq;
%% ===========================================================================

% Loop over all defined cable loops
for l=1:length(loops),
fprintf('\n***** Loop #%d ***** \n',l);

% Loop over all PSD sets
for psdbase=psdset,
ETSImasks.cab = sprintf('%s%sM%d', mask.base, mask.ext.cab, psdbase);
ETSImasks.ex = sprintf('%s%sM%d', mask.base, mask.ext.ex, psdbase)
switch psdbase
    case 1,
        capval=-51.44; % reduce power by cutting the top PSD values 
        tmp_tfplan.PSD.HAM.active = 0; % Our HAM band  stuff do not work
                                       % correctly with PBO
    case 2,
        capval=-55.5; 
        tmp_tfplan.PSD.HAM.active = 0;
    otherwise error('Unknown psdbase');
end;
    
% Loop over all Noises
for a=1:length(FSANnoise),
    switch FSANnoise{a}
        case {'0'},
            noisetype='Calculate'; 
            ETSImask = ETSImasks.ex;
            freqplan = [freqplanbase '-Ex'];
	case {'A','B','C'},
            noisetype=sprintf('Model %s (FTTCab VDSL)',FSANnoise{a});
            ETSImask = ETSImasks.cab;
            freqplan = [freqplanbase '-Cab'];
	case {'D','E','F'},
            noisetype=sprintf('Model %s (FTTEx VDSL)',FSANnoise{a});	
            ETSImask = ETSImasks.ex;
            freqplan = [freqplanbase '-Ex'];
	otherwise  error('Unknown noise model');
    end;
    ex.param.FSANNoiseModel=noisetype;   

    % Now the parameters for the experiment is set let uss get what the user wished for
    masktf = getList(ex.tfplist,ETSImask);

    % Extract the plan
    downstream_val=str2num(masktf.PSD.downstream(9:findstr(masktf.PSD.downstream,']')));
    upstream_val  =str2num(masktf.PSD.upstream(9:findstr(masktf.PSD.upstream,']')));

    % Definitions of VDSL PSDs (frequency plans)
    switch freqplan
        case { 'NBand-Cab', 'NBand-Ex'}, 
            fmin=fq(1); fmax=fq(end);
            fqstr=sprintf('%g ',fq/1e6);
            duplexname=sprintf('%d-band-plan %sMHz (ETSI mask %s)',...
                length(fq)-1,fqstr, ETSImask);
            UFplan=makeUFplan(fq,fqdir(1:length(fq)-1));
            [tmp_upstream,tmp_downstream]= setupPSDplan(UFplan,...
                upstream_val, downstream_val, 'Log-Linear', -120,capval);
        otherwise
            estr = sprintf('Unknown freqency plan: %s',freqplan); error(estr);
    end; %switch

    % Set up the new PSD mask
    tmp_tfplan.PSD.active.upstream=[fmin fmax];
    tmp_tfplan.PSD.active.downstream=[fmin fmax];
    tmp_tfplan.PSD.downstream=['calcPSD([', num2str(tmp_downstream),'],''Log-Linear'',ex.param.frequency.f)'];
    tmp_tfplan.PSD.upstream  =['calcPSD([', num2str(tmp_upstream),'],''Log-Linear'',ex.param.frequency.f)'];
    ex.tfplist=setList(ex.tfplist,tmp_tfplan.name,tmp_tfplan);

    % plot PSD mask
    if l==loops(1),
        switch FSANnoise{a}
            case {'A','E'}, 
                if doplots,
                    fax.min=3e3;
                    fax.max=fmax+1e6;
                    figure;
                    axes('Units','pixels');
                    plotTFplan(tmp_tfplan,'Linear',fax);
                    drawnow;
                end; %doplots
                fprintf('HAMband flag = %.0f (%s)\n',...
                    tmp_tfplan.PSD.HAM.active,ex.param.HAMBandName);
                fprintf('Cap val = %5.2f\n',capval);
        end;
    end;

    % Print experiment setup
    fprintf('\n*****\n')
    fprintf('%s\n',duplexname);
    fprintf('Uses 4x5 dist, AWGN -140 and Noise %s\n',noisetype)

    fprintf('Goal bitrates:\n');
    fprintf('%6.2f\t', Grates(distorder,2)/1e6) % Downstream
    fprintf('\n')
    fprintf('%6.2f\t', Grates(distorder,1)/1e6) % Upstream
    fprintf('\n')
    fprintf('Modem distances:\n')
    fprintf('%6.0f\t', Gdists(distorder))
    fprintf('\n')
    
    % Set up topology
    nodenames=Gnames(distorder,:);
    dd= [dists(1) diff(dists)];
    switch mss
        case 1,
            tt.name=sprintf('ETSI Mixed services from Cabinet');
            switch l
                case 1, cable='BT_dwug';  % TP100
                    tt.topology=[{0  '' 'CO' ''}; {200  '' 'C' 'TP100'}];
                    for q=1:length(dists),
                        tt.topology=[tt.topology; ...
                                {dd(q)  cable sprintf('%s',nodenames(q,:)) sprintf('%.0f m',dd(q))}; ];
                    end;
                    tt.traffic=[{}];
                    for q=1:length(dists),
                        tt.traffic=[tt.traffic; ...
                                {2 q+2 'VDSL' nomodemspernode}; ];
                    end;
                otherwise  estr=sprintf('Loop #d not implemented\n',l); error(estr);
            end;
        case 2,
            tt.name=sprintf('ETSI Mixed services from Exchange');
            switch l
                case 1, cable='BT_dwug';  % TP100
                    tt.topology=[{0  '' 'CO' ''}];
                    for q=1:length(dists),
                        tt.topology=[tt.topology; ...
                                {dd(q)  cable sprintf('%s',nodenames(q,:)) sprintf('%.0f m TP100',dd(q))}; ];
                    end;
                    tt.traffic=[{}];
                    for q=1:length(dists),
                        tt.traffic=[tt.traffic; ...
                                {1 q+1 'VDSL' nomodemspernode}; ];
                    end;
                case 2, cable='KPN_L1';  % TP150
                    tt.topology=[{0  '' 'CO' ''}];
                    for q=1:length(dists),
                        tt.topology=[tt.topology; ...
                                {dd(q)  cable sprintf('%s',nodenames(q,:)) sprintf('%.0f m TP150',dd(q))}; ];
                    end;
                    tt.traffic=[{}];
                    for q=1:length(dists),
                        tt.traffic=[tt.traffic; ...
                                {1 q+1 'VDSL' nomodemspernode}; ];
                    end;
                otherwise  estr=sprintf('Loop #d not implemented\n',l); error(estr);
            end;

    end;
    ex.tt=tt;				% Define the experiment tt structure

    for p=1:length(pbomethods),
        tmp_tfplan.PSD.PBO.method	= pbomethods{p};
        tmp_tfplan.PSD.PBO.param.len	= pboparam.len{p};	
        tmp_tfplan.PSD.PBO.param.freq	= pboparam.freq{p};	
        ex.tfplist=setList(ex.tfplist,tmp_tfplan.name,tmp_tfplan);

        % Print what PBO we are using
        fprintf('\n-----\n');
        fprintf('PBO method: %s\nPBO parameters f =',pbomethods{p});
        fprintf(' %g',pboparam.freq{p});
        fprintf(', l =');
        fprintf(' %g',pboparam.len{p});
        fprintf(' (maxl = %g)\n',pboparam.maxlen);

        
        % Evaluate this experiment
        result = evalExperiment;
        [bitrate_LT, bitrate_NT, dummy, dummy2]=calcFSANresult(ex,result);

        % Save the experiment for later use
        simures(p).ex=ex;
        simures(p).result=result;
        simures(p).LT=bitrate_LT;
        simures(p).NT=bitrate_NT;
        simures(p).model=ex.param.FSANNoiseModel;
        simures(p).scenario=tt.name;

        symrate = min(bitrate_NT,bitrate_LT);
        tab(q,:)=[ bitrate_LT, bitrate_NT];

        if printpower,
            f=ex.param.frequency.f;
            fi=(1:30000).*1e3;
            fprintf ('Actual total Power:\nNT [dBm]:');
            dl=0;
            for rno=1:length(result),
                PSD_NT=result(rno).NT.Tx_signal;
                PSD_LT=result(rno).LT.Tx_signal;
                PSD_NTi=interp1(f,PSD_NT,fi);
                fprintf (' %2.2f',...
                    10*log10(sum(PSD_NTi))+10*log10(fi(4)-fi(3)));
                dl=dl+ ex.tt.topology{result(rno).Modem.NT_Node,1};
                fprintf(' @%2.0fm,',dl);
            end;
            PSD_LTi=interp1(f,PSD_LT,fi);
            fprintf ('\nLT=%2.2f dBm\n',...
                10*log10(sum(PSD_LTi))+10*log10(fi(4)-fi(3)));
        end;
        
        % Print experiment results
        fprintf('\nBitrates:\n')
        fprintf('%6.2f\t', bitrate_NT)
        fprintf ('\n');
        fprintf('%6.2f\t', bitrate_LT)
        fprintf ('\n');
        diffd = bitrate_NT' ./ (Grates(distorder,2)./1e6 );
        diffu = bitrate_LT' ./ (Grates(distorder,1)./1e6 );
        diffa = min(diffd,diffu); % actual diff
        fprintf('Percentage of goal:\n')
        fprintf('%6.2f\t', diffa*100)
        fprintf ('\n');
        fprintf('(Sum bitrate downstrean = %3.2f, upstrean = %3.2f; Total = %3.2f)\n', sum(bitrate_NT), sum(bitrate_LT), sum([bitrate_NT bitrate_LT]))
        
        % Plot experiment results
        if doplots,
            figure;
            txno=length(simures(p).result);
            txval=[];
            dl=[];
            for txi=1:txno,
                txval(txi,:)=10*log10(simures(p).result(txi).NT.Tx_signal);
                dl=[dl ex.tt.topology{simures(p).result(txi).Modem.NT_Node,1}];
            end;
            plot(f,txval);
            PBOinfo=simures(p).ex.tfplist(end).PSD.PBO;

            tstr=[sprintf('PBO method %s; l=', PBOinfo.method), ...
                    sprintf('%d ', PBOinfo.param.len) 'f=', ...
                    sprintf('%g ', PBOinfo.param.freq./1e6), ...
                    sprintf('(maxl = %g)',pboparam.maxlen) ];
            title(tstr);
            legend(num2str(cumsum(dl')))
            drawnow;
        end;
        
    end; % PBO methods

end; % FSANnoise

end; % psd set

end; %loops