function [PSD_up] = calcUPBO(PSD_up,PBO,Att_mtr,start_node,stop_node,Ref_Topology,Ref_Noise);
%% ===========================================================================
% calcUPBO -     Calculates upstream PBO
%
%        INPUT
%        -----
%        PSD_up    :    PSD template (upstream) (Tx signal)
%        PBO       :    PBO struct (method  and param)
%        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; V2.0.5 (2002-11)
%% ===========================================================================

%% ===========================================================================
% Copyright (C):                                        
%       1998-1999 by Telia Research AB, Lulea, Sweden;                
%       2000-2003 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,v 3.5 2003/11/03 12:33:22 tono Exp $
%% ===========================================================================
% 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
%% ===========================================================================
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
expNoise=ex.param.expNoise; 
pbo_index_up=find(PSD_up>10.^expNoise);

switch PBO.method

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

 case 'None', 	
  % No PBO method
  % Do Nothing (i.e., input PSD_up is the same as output PSD_up)!
  Zterm=Zterm;

%% ===========================================================================
 case 'RefRef',	
  % Reference PBO method [7]
  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;
  LOSS = kl*sqrt(f/1e6);
  a=PBO.param.refa;
  b=PBO.param.refb;

  TxPSD=-140*ones(1,length(f)); %10*log10(PSD_up);
  
  nobands=length(PBO.param.upf_start);
  if length(a) == 1; % If we have only one parameter for all bands
      PSDref = a+b*sqrt(f);
      TxPSD  = PSDref + LOSS;
  else
      for band = 1:nobands;
          fix=find(f>PBO.param.upf_start(band) & f<PBO.param.upf_stop(band));
          PSDref = a(band)+b(band)*sqrt(f(fix)/1e6);
          TxPSD(fix)=PSDref + LOSS(fix);
      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}=PBO.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=PBO.param.plan.up.fc+PBO.param.plan.up.bw/2;
  low =PBO.param.plan.up.fc-PBO.param.plan.up.bw/2;
  E_up=ones(1,length(f));
  for band=1:length(PBO.param.len),
     Ref_Topology{2}{1}=PBO.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=PBO.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=PBO.param.plan.up.fc+PBO.param.plan.up.bw/2;
  low =PBO.param.plan.up.fc-PBO.param.plan.up.bw/2;
  E_up=ones(1,length(f));
  
  for band=1:length(PBO.param.len),
      l_ref=PBO.param.len(band);
      f_ref=PBO.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=PBO.param.freq(1); % find the constant
  f_index=max(find(f<=f_ref));
  Ref_Topology{2}{1}=PBO.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=PBO.param.plan.up.fc+PBO.param.plan.up.bw/2;
  low =PBO.param.plan.up.fc-PBO.param.plan.up.bw/2;
  E_up=ones(1,length(f));
  Ref_Topology{2}{1}=PBO.param.maxlen;
  Ref_Att=getIL(ex,1,2,Ref_Topology,f,Zterm);
  for band=1:length(PBO.param.freq),
      f_ref=PBO.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 ready yet!
  % Flat spectrum constant average log method (ITU Q4/15 D.815, Geneva 2000)
  Tx=10*log10(PSD_up);
  Rx=-60*ones(1,length(f)); % ToDo: Get correct values for Rx!
  high=PBO.param.plan.up.fc+PBO.param.plan.up.bw/2;
  low =PBO.param.plan.up.fc-PBO.param.plan.up.bw/2;
  E_up=ones(1,length(f));
  for band=1:length(PBO.param.k),
      k_ref=PBO.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=PBO.param.plan.up.fc+PBO.param.plan.up.bw/2;
  low =PBO.param.plan.up.fc-PBO.param.plan.up.bw/2;
  E_up=ones(1,length(f));
  len_att=sumLen(start_node,stop_node)./1000;

  for band=1:length(PBO.param.freq),
      q_ref=PBO.param.q(band);    % get the boosting parameter
      f_ref=PBO.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',PBO.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)));