function [Rate_LT, Rate_NT, Margin_LT, Margin_NT] = calcResultDMT(tfplan,res,lc,f);
%% ===========================================================================
%calcDMTRate - Calculates the total rate for a DMT system
% Parameter:    tfplan   	The TF plan used for calculating the bit rate
% Parameter:    res	     	Result structure
% Parameter:    lc	        The linecode
% Parameter:    f	        Frequency axis
%
% Returns:      Rate__LT    	Resulting bitrates at LT side (in Mbit)
% Returns:      Rate__NT    	Resulting bitrates at NT side (in Mbit)
% Returns:      Margin_LT 
% Returns:      Margin_NT
% Example(s):
%
% Algorithmic details:
%	The calculation of the rate is basically an approximative
%	calculation of an integral, with b and f as the 
%	y- and x-axis. The algorithm approximates the integral by 
%	taking the average value between two points multiplied 
%  	with the frequency difference of the two points. 
%
% Reference:
%    Driton Statovci, Rickard Nilsson, and Tomas Nordström, 
%    "Generic detection model for DMT based modems", Temporary
%    Document TD 23r2, 034t23r2, ETSI STC TM6, Sophia Antipolis, 
%    France, 17-21, Nov, 2003
%% ===========================================================================

%% ===========================================================================
% Copyright (C):                                        
%       1999-2000 by Telia Research AB, Lulea, Sweden;                
%       2001-2003 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                         All rights reserved.
% Project       : FTW's xDSLsimu
% Author(s)     : Daniel Bengtsson (Daniel.J.Bengtsson@Telia.se)
%               : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%
% CVS:       $Id: calcResultDMT.m,v 3.7 2003/12/23 10:35:07 tono Exp $
%% ===========================================================================
% Change History
%      1999-08-13 (DaB)  Created
%      1999-10-11 (ToNo) Reimplemented, now based on the new calcTheoRate
%      1999-11-22 (DaB)  Margin output added
%      2000-01-12 (DaB)  Using dmt function as bitloading function
%      2000-04-03 (DaB)  CalcRate removed
%      2000-04-04 (DaB)  Signal_margin added
%      2000-04-06 (UvAn) Passing ftplan, lc and f into the function; 
%                        removed global ex
%      2000-04-11 (DaB)  Added sync switch
%      2001-06-20 (Bosse) Octave port (first time)
%      2001-07-03 (ToNo) Calculate bit-rates for active tones
%      2001-09-12 (ToNo) Add support for "Maximum aggregate transmit power"
%      2001-12-14 (ToNo) Cleaning up the code (remove unused things)
%      2002-08-26 (ToNo) Made it possible to change water-filling method
%      2002-10-02 (ToNo) Revised the input parameters to the water-fill alg.'s
%      2003-03-20 (ToNo) Adjusted for template vs. peak mask input to WF
%      2003-12-02 (ToNo) Added support for lc dependent modem noise
%      2003-12-16 (ToNo) Split background noise into bgNoise and modemNoise
%      2003-12-18 (DrSt) Scale downstream and upstream powers depending on
%                        cyclic extension
%      2003-12-18 (DrSt) Renamed BitLoading (BL) algorithms
%% ===========================================================================

if nargin<4,
    error('Not enough input arguments.');
end

% Generic parameters defined in lc list
efficiencyLoss	= lc.param.efficiencyLoss;
SNRmax	        = lc.param.SNRMax;
codingGain      = lc.param.codingGain;
SNRloss	        = lc.param.SNRloss; 
refSNR          = lc.param.refSNR;
xtalk_margin    = lc.param.xtalk_margin;
signal_margin   = lc.param.signal_margin;
PxUp            = lc.param.PxUp;
PxDown          = lc.param.PxDown;

% DMT parameters
df              = lc.param.dmt.df;
srate           = lc.param.dmt.srate;
constellation   = lc.param.constellation;
ds              = lc.param.dmt.ds;        % Get carriers, tones, cp, cs
us              = lc.param.dmt.us;


newf= df:df:df*ds.carriers; % The frequency of the DMT tones

% Get some time and frequency parameters
TimeDivision_up         = tfplan.timeDivision.up;
TimeDivision_down       = tfplan.timeDivision.down;
sync                    = tfplan.timeDivision.sync;
f_start_NT              = tfplan.PSD.active.downstream(1);
f_stop_NT               = tfplan.PSD.active.downstream(2);
f_start_LT              = tfplan.PSD.active.upstream(1);
f_stop_LT               = tfplan.PSD.active.upstream(2);

% Pick up results from the result structure
Rx_NT=res.NT.Rx_signal;
Tx_NT=res.NT.Tx_signal;
Rx_LT=res.LT.Rx_signal;
Tx_LT=res.LT.Tx_signal;
if (sync==1),
   Tot_NT=res.NT.Tot_noise.down;
   Tot_LT=res.LT.Tot_noise.up;
else
   Tot_NT=res.NT.Tot_noise.up;
   Tot_LT=res.LT.Tot_noise.down;
end;

% The total noise is increased by xtalk_margin dB:s
% (which is not in the ShannonGap).
noise_LT=Tot_LT.*10.^(xtalk_margin/10);
noise_NT=Tot_NT.*10.^(xtalk_margin/10);

% Add modem noise
if exist('lc.param.modemNoise.LT')~=1
    noise_LT = noise_LT + 10.^(eval([lc.param.modemNoise.LT ';'])/10);
    noise_NT = noise_NT + 10.^(eval([lc.param.modemNoise.NT ';'])/10);
end;

% Prepare for water-filling
H_LT=interpolate(f,Rx_LT./Tx_NT,newf);
N_LT=interpolate(f,noise_LT,newf);
H_NT=interpolate(f,Rx_NT./Tx_LT,newf);
N_NT=interpolate(f,noise_NT,newf);
PSD_NT=interpolate(f,Tx_NT,newf);
PSD_LT=interpolate(f,Tx_LT,newf);

if lc.name(1) == 'A'
    % Limit the total power (this is a formula for ADSL)
    mDown=ceil(10*(-3.65+10*log10(length(ds.tones)+1)))/10;
    mDown=min(PxDown,mDown);
    mDown=10.^(mDown./10)./df;
    mUp=  ceil(10*(-1.65+10*log10(length(us.tones)+1)))/10;
    mUp=min(PxUp,mUp);
    mUp=10.^(mUp./10)./df;
else
    mDown=10.^(PxDown./10)./df; 
    mUp  =10.^(PxUp./10)./df;
end;

% Note that the channel magnitude response here is based on the
% insertion loss, as the PSD is the mask at the input of the line.

% Reduce the power according to the power lost to due the cyclic extension.
if ~isfield(lc.param,'skippowerfix')
    us.ce=us.cp + us.cs;
    ds.ce=ds.cp + ds.cs;
    mUp=mUp*((2*us.carriers)/(2*us.carriers+us.ce));
    mDown=mDown*((2*ds.carriers)/(2*ds.carriers+ds.ce));
end

% Do DMT water-filling
switch tfplan.fixBitrate.name;
 case {'None','RBL'} % RBL-Rounded Bit-Loading
  % The original rate adaptive water-filling, which does rounding and
  % can therefore output to much power!
  Gap=10.^((refSNR+SNRloss+signal_margin-codingGain)/10);
  [b_LT, E_LT, g_LT, M_LT] = bitloadingRounded(sqrt(abs(H_LT)),N_LT,mUp,Gap,PSD_NT,constellation,us.tones);
  [b_NT, E_NT, g_NT, M_NT] = bitloadingRounded(sqrt(abs(H_NT)),N_NT,mDown,Gap,PSD_LT,constellation,ds.tones);
case {'GABL_RA'}
  %Levin-Campello RA implemenation of Gain Adjusted Bit-Loading(GABL)
  Gap=10.^((refSNR+SNRloss+signal_margin-codingGain)/10);
  LC_type=0;        % Levin-Campello RA (rate adaptive)
  bitrate_up=0;
  bitrate_down=0;
  [b_LT, E_LT, g_LT, M_LT] = bitloadingGA(sqrt(abs(H_LT)),N_LT,mUp,Gap,...
     PSD_NT,constellation,us.tones,LC_type,bitrate_up, srate, res.Modem);
  [b_NT, E_NT, g_NT, M_NT] = bitloadingGA(sqrt(abs(H_LT)),N_NT,mDown,Gap,...
     PSD_LT,constellation,ds.tones,LC_type,bitrate_down, srate, res.Modem);
  
 case {'GABL_MA'} 
  %Levin-Campello MA implemenation of Gain Adjusted Bit-Loading(GABL)
  Gap=10.^((refSNR+SNRloss+signal_margin-codingGain)/10);
  LC_type=1;        % Levin-Campello MA (margin adaptive)
  bitrate_up=tfplan.fixBitrate.param.target_us;
  bitrate_down=tfplan.fixBitrate.param.target_ds;
  [b_LT, E_LT, g_LT, M_LT] = bitloadingGA(sqrt(abs(H_LT)),N_LT,mUp,Gap,...
     PSD_NT,constellation,us.tones,LC_type,bitrate_up, srate, res.Modem);
  [b_NT, E_NT, g_NT, M_NT] = bitloadingGA(sqrt(abs(H_LT)),N_NT,mDown,Gap,...
     PSD_LT,constellation,ds.tones,LC_type,bitrate_down, srate, res.Modem);
 
 otherwise
  estr = sprintf('Error in calcResultDMT, unknown WF method: %s\n',...
                 tfplan.fixBitrate.name);
  error(estr);
end

%Base bitrate on the active tones
Rate_LT=sum(b_LT(us.tones)).*srate./1e6.*TimeDivision_up.*(1-efficiencyLoss);
Rate_NT=sum(b_NT(ds.tones)).*srate./1e6.*TimeDivision_down.*(1-efficiencyLoss);

Margin_LT=mean(M_LT)+xtalk_margin+signal_margin;
Margin_NT=mean(M_NT)+xtalk_margin+signal_margin;