function [PSD_up] = calcUPBO(PSD_up,UPBO,Att_mtr,start_node,stop_node,Ref_Topology,Ref_Noise)
%% ===========================================================================
% calcUPBO -     Calculates upstream PBO
%
%        INPUT
%        -----
%        PSD_up    :    PSD template (upstream) (Tx signal)
%        UPBO      :    UPBO struct (method, param, and band)
%        Att_mtr   :    Matrix with Attenuation in dB/km
%        start_node:    Start node for PBO modem
%        stop_node :    Stop node for PBO modem
%        Ref_Topology : Reference Topology for PBO calculation
%        Ref_Noise :    Reference Noise for PBO method
%
%        OUTPUT
%        ------
%        PSD_up    :    Transmit PSD after PBO
%
% Algorithmic details:
%
% Reference:
%     [1] ETSI Document TS 101 270-2, Transmission and Multiplexing
%         (TM); Access transmission systems on metal access cables; Very
%         high-speed Digital Subscriber Line (VDSL); Part 2: Tranceiver
%         Specification; V0.1.2 (2000-03)
%     [2] FTW's xDSLsimu manual
%     [3] FSAN VDSL working group, "Power-backoff methods for VDSL",
%         ETSI TM6 983T17A0, Lulea, Sweden June 1998
%     [4] ITU Q4/15 FI-102, Fiji 2000
%     [5] ITU Q4/15 FI-114, Fiji 2000
%     [6] ITU Q4/15 D.815, Geneva 2000
%     [7] ETSI Document TS 101 270-1, Transmission and Multiplexing
%         (TM); Access transmission systems on metal access cables; Very
%         high-speed Digital Subscriber Line (VDSL); Part 1: Functional
%         Requirements; V1.3.1 (2003-07)
%% ===========================================================================

%% ===========================================================================
% Copyright (C):                                        
%       1998-1999 by Telia Research AB, Lulea, Sweden;                
%       2000-2009 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                         All rights reserved.
% Project       : FTW's xDSLsimu
% Author(s)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%               : Daniel Bengtsson (Daniel.J.Bengtsson@Telia.se)
%               : Bo Engstrom (bosse@upzide.com)
%
% CVS:       $Id: calcUPBO.m 752 2009-01-02 13:03:52Z tono $
%% ===========================================================================
% Change History
%      1998-12-16 (DaB)  Created     
%      1998-12-22 (ToNo) Added and updated comments
%      1999-02-04 (DaB)  Modified to allow multiple cables 
%      1999-03-10 (DaB)  Added new PBO methods
%      1999-03-11 (DaB)  New input arguments
%      1999-08-11 (DaB)  Rewritten for new structure
%      1999-09-30 (DaB)  Rewritten for new ex struct
%      1999-11-16 (ToNo) Implementing the multiple reference length method
%      2000-03-23 (ToNo) Implementing the multiple reference freq method
%      2000-04-04 (ToNo) Implementing a number of multiple constant PBO methods
%      2001-05-30 (ToNo) Fixed a bug in how to calculate the
%                        reference attenuation for various lengths
%      2001-06-18 (Bosse) Octave port
%      2002-07-01 (PeKa) Octave to Matlab compatibility port
%      2003-01-02 (ToNo) Added Reference PBO method
%      2003-01-02 (ToNo) Changed parameters for evalPBO
%                        Instead of a scaling we now use the actual PSD
%      2003-01-03 (ToNo) Split upstream and downstream PBO
%      2003-01-07 (ToNo) Fixed that octave do not want empty case in switch
%      2003-04-11 (ToNo) Fixed bug in RefRef (loss should be positive)
%      2003-04-22 (ToNo) Added the limit on kl (>= 1.8) for RefRef
%      2003-11-03 (ToNo) Lists now use cell arrays for both Octave and Matlab
%      2005-03-29 (ToNo) Consolidated the way to specify the upstream bands
%      2005-04-02 (ToNo) Fixed a bug in single band definition of RefRef
%      2005-04-12 (ToNo) Added an ideal RefRef method
%      2005-05-10 (ToNo) Changed how we initialise the result in RefRef*
%      2008-10-08 (ToNo) Allowed turing of RefRef UPBO by setting
%                        parameter UPBO.param.refa to zero
%      2008-10-08 (ToNo) Added support for "final el" (Mgmt 'el' override)
%% ===========================================================================
global ex;

Zterm=ex.param.Zterm;            
f=ex.param.frequency.f; 	% Frequency axis

% Determine the minimum Tx value to be used
minValue_up=min(PSD_up);

% Determine where there are any Tx at all
% FIXME this is just a bad hack for old PBO methods,
% we should use UPBO.band instead
pbo_index_up=find(PSD_up>10.^-12);

switch UPBO.method

%% ===========================================================================

 case {'None', 'none'}
  % No PBO method
  % Do Nothing (i.e., input PSD_up is the same as output PSD_up)!
  Zterm=Zterm; % A no-operation to reduce octave complaints

%% ===========================================================================
 case 'RefRefIdeal',	
  % Reference PBO method using ideal loss knowledge and no 1.8 limit
  LOSS =-getIL(ex,start_node,stop_node, ex.tt.topology, f, ex.param.Zterm);
 
  a=UPBO.param.refa;
  b=UPBO.param.refb;

  TxPSD=PSD_up; % Initialize the output, esp. in the out-of-band part

  nobands=length(UPBO.band(:,1));
  if length(a) == 1; % If we have only one parameter for all bands
      PSDref = a+b*sqrt(f/1e6);
      TxPSD  = PSDref + LOSS;
  else
      for band = 1:nobands;
          if a(band) ~= 0,
              fix=find(f>UPBO.band(band,1) & f<UPBO.band(band,2));
              PSDref = a(band)+b(band)*sqrt(f(fix)/1e6);
              TxPSD(fix)=PSDref + LOSS(fix);
          end;
      end;
  end;
  PSD_up = min(PSD_up,10.^(TxPSD/10)); % Do not get over the mask
   
%% ===========================================================================
 case 'RefRef',	
  % Reference PBO method [7]
  % Estimate an electrical length kl and model LOSS according to sqrt(f)
  freq = 1e6:1e6:12e6;
  loss = -getIL(ex,start_node,stop_node,ex.tt.topology,freq,ex.param.Zterm);
  kl = min(loss./sqrt(freq/1e6));
  if kl<1.8, kl=1.8; end;
  if isfield(UPBO.param,'finalel'), kl = UPBO.param.finalel; end;
  LOSS = kl*sqrt(f/1e6);
  a=UPBO.param.refa;
  b=UPBO.param.refb;

  TxPSD=PSD_up; % Initialize the output, esp. in the out-of-band part
  
  nobands=length(UPBO.band(:,1));
  if length(a) == 1; % Only one parameter for all bands?
      PSDref = a+b*sqrt(f/1e6);
      TxPSD  = PSDref + LOSS;
  else
      for band = 1:nobands;
          if a(band) ~= 0,
              fix=find(f>UPBO.band(band,1) & f<UPBO.band(band,2));
              PSDref = a(band)+b(band)*sqrt(f(fix)/1e6);
              TxPSD(fix)=PSDref + LOSS(fix);
          end;
      end;
  end;
  PSD_up = min(PSD_up,10.^(TxPSD/10)); % Do not get over the mask
   
%% ===========================================================================
 
 case 'RefLen',	
  % Reference length method
  att=Att_mtr{start_node,stop_node};
  Ref_Topology{2}{1}=UPBO.param.len(1);
  Ref_Att=getIL(ex,1,2,Ref_Topology,f,Zterm);
  tmp=Ref_Att-att;
  tmp=min(tmp,zeros(size(f)));
  E_up=ones(1,length(f)).*10.^(tmp/10);
  PSD_up(pbo_index_up)=PSD_up(pbo_index_up) .* E_up(pbo_index_up);
  
 case 'RefLenM',	
  % Reference length method with multiple references
  att=Att_mtr{start_node,stop_node};
  high=UPBO.band(:,2);
  low =UPBO.band(:,1);
  E_up=ones(1,length(f));
  for band=1:length(UPBO.param.len),
     Ref_Topology{2}{1}=UPBO.param.len(band);
     Ref_Att=getIL(ex,1,2,Ref_Topology,f,Zterm);
     tmp=Ref_Att-att;
     tmp=min(tmp,zeros(size(f)));
     ix=find(f<high(band) & f>low(band));
     E_up(ix)=10.^(tmp(ix)/10);
  end;
  PSD_up(pbo_index_up)=PSD_up(pbo_index_up) .* E_up(pbo_index_up);
   
%% ===========================================================================

 case 'RefFEXT',	
  % Reference FEXT method
  l_ref=UPBO.param.len(1);
  att=Att_mtr{start_node,stop_node};
  len_att=sumLen(start_node,stop_node);
  Ref_Topology{2}{1}=l_ref;
  Ref_Att=getIL(ex,1,2,Ref_Topology,f,Zterm);
  tmp=Ref_Att+10*log10(l_ref) ...
      -att-10*log10(len_att);
  tmp=min(tmp,zeros(size(f)));
  E_up=ones(1,length(f)).*10.^(tmp/10);
  PSD_up(pbo_index_up)=PSD_up(pbo_index_up) .* E_up(pbo_index_up);

 case 'ConstantRefFEXTM'	
  % Flat spectrum reference FEXT (ITU FI-114)
  att=Att_mtr{start_node,stop_node};
  len_att=sumLen(start_node,stop_node);
  high=UPBO.band(:,2);
  low =UPBO.band(:,1);
  E_up=ones(1,length(f));
  
  for band=1:length(UPBO.param.len),
      l_ref=UPBO.param.len(band);
      f_ref=UPBO.param.freq(band); % find the level at the reference frequency 
      f_index=max(find(f<=f_ref));
      Ref_Topology{2}{1}=l_ref;
      Ref_Att=getIL(ex,1,2,Ref_Topology,f,Zterm);
      tmp=Ref_Att(f_index)+10*log10(l_ref) ...
	  -att(f_index)-10*log10(len_att);
      tmp=min(tmp,zeros(size(f)));
      ix=find(f<high(band) & f>low(band));
      E_up(ix)=10.^(tmp(ix)/10);
  end;
  PSD_up(pbo_index_up)=PSD_up(pbo_index_up) .* E_up(pbo_index_up);

%% ===========================================================================

 case {'RefFreq' , 'Constant'} 
  % Reference frequency method (Constant power)
  att=Att_mtr{start_node,stop_node};
  f_ref=UPBO.param.freq(1); % find the constant
  f_index=max(find(f<=f_ref));
  Ref_Topology{2}{1}=UPBO.param.maxlen;
  Ref_Att=getIL(ex,1,2,Ref_Topology,f,Zterm);
  tmp=Ref_Att(f_index)-att(f_index);
  tmp=min(tmp,zeros(size(f)));
  E_up=ones(1,length(f)).*(10.^(tmp/10));
  PSD_up(pbo_index_up)=PSD_up(pbo_index_up) .* E_up(pbo_index_up);

 case {'RefFreqM', 'ConstantM'}
  % Reference frequency (Constant) method with multiple references
  att=Att_mtr{start_node,stop_node};
  high=UPBO.band(:,2);
  low =UPBO.band(:,1);
  E_up=ones(1,length(f));
  Ref_Topology{2}{1}=UPBO.param.maxlen;
  Ref_Att=getIL(ex,1,2,Ref_Topology,f,Zterm);
  for band=1:length(UPBO.param.freq),
      f_ref=UPBO.param.freq(band); % find the constant
      f_index=max(find(f<=f_ref));
      tmp=Ref_Att(f_index)-att(f_index);
      tmp=min(tmp,zeros(size(f)));
      ix=find(f<high(band) & f>low(band));
      E_up(ix)=10.^(tmp(ix)/10);
  end;
  PSD_up(pbo_index_up)=PSD_up(pbo_index_up) .* E_up(pbo_index_up);

%% ===========================================================================

 case 'RefNoise'	
  % Reference Noise method	
  att=Att_mtr{start_node,stop_node};
  len_att=sumLen(start_node,stop_node)./1000;
  tmp=10*log10(Ref_Noise)-att-20*log10(f./1e6)-10*log10(len_att)-...
      ex.param.XTlevel.FEXT;
  tmp=min(tmp,zeros(size(f)));
  E_up=ones(1,length(f)).*(10.^(tmp/10));
  PSD_up(pbo_index_up)=PSD_up(pbo_index_up) .* E_up(pbo_index_up);

%% ===========================================================================
  
 case 'ConstantAvgLogM'	% NOTE/WARNING - This is not tested yet!
  % Flat spectrum constant average log method (ITU Q4/15 D.815, Geneva 2000)
  Tx=10*log10(PSD_up);
  att=Att_mtr{LT_node,NT_node};
  Rx=Tx+att;
  Rx=(10.^(Rx/10));
  high=UPBO.band(:,2);
  low =UPBO.band(:,1);
  E_up=ones(1,length(f));
  for band=1:length(UPBO.param.k),
      k_ref=UPBO.param.k(band);
      ix=find(f<high(band) & f>low(band));
      k=trapz(f(ix),log2(Rx(ix)));
      level=2^((k_ref-k)/(high(band)-low(band))); %%% Should this be 10*log10 ?
      E_up(ix)=level*Tx(ix);
  end;
  PSD_up(pbo_index_up)=PSD_up(pbo_index_up) .* E_up(pbo_index_up);
  
 case 'ConstantRefNoiseM'	
  % Flat spectrum reference noise (ITU Q4/15 FI-114, Fiji 2000)
  att=Att_mtr{start_node,stop_node};
  high=UPBO.band(:,2);
  low =UPBO.band(:,1);
  E_up=ones(1,length(f));
  len_att=sumLen(start_node,stop_node)./1000;

  for band=1:length(UPBO.param.freq),
      q_ref=UPBO.param.q(band);    % get the boosting parameter
      f_ref=UPBO.param.freq(band); % find the level at the reference frequency 
      f_index=max(find(f<=f_ref));
      tmp=10*log10(q_ref)+10*log10(Ref_Noise) ...
	  -att(f_index)-20*log10(f_ref./1e6)-10*log10(len_att)-...
          ex.param.XTlevel.FEXT;
      tmp=min(tmp,zeros(size(f)));
      ix=find(f<high(band) & f>low(band));
      E_up(ix)=10.^(tmp(ix)/10);
  end;
  PSD_up(pbo_index_up)=PSD_up(pbo_index_up) .* E_up(pbo_index_up);

%% ===========================================================================

 otherwise
  sprintf('Unknown UPBO method: %s',UPBO.method)
end; %switch


% Set min Tx to the min of the original Tx
PSD_up(find(PSD_up<minValue_up))=minValue_up*...
    ones(size(find(PSD_up<minValue_up)));