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:
%    FSAN xDSL simulation tool manual
%% ===========================================================================

%% ===========================================================================
% Copyright (C) 1999 by Telia Research AB, Lulea, Sweden; All rights reserved. 
% Project       : FSAN duplex model
% Author(s)     : Daniel Bengtsson (Daniel.J.Bengtsson@Telia.se)
%               : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%
% CVS:       $Id: calcResultDMT.m,v 1.4 2000/05/03 11:29:31 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
%% ===========================================================================

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

% 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;
constellation   = lc.param.constellation;
xtalk_margin    = lc.param.xtalk_margin;
signal_margin   = lc.param.signal_margin;
Px              = lc.param.Px;

df              = lc.param.dmt.df;
carriers        = lc.param.dmt.carriers;

newf= df:df:df*carriers;


reqSNR = 10*log10(10.^((refSNR+SNRloss-codingGain+signal_margin)/10).*(2.^(constellation)-1)); 

% Get sometime 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);

% Calculate SNR
Rx_NT=res.NT.Rx_signal;
Tx_NT=res.NT.Tx_signal;
tmp=find(Tx_NT>10.^-14);
LT_mask=10.^14.*ones(size(f));
LT_mask(tmp)=ones(size(tmp));
Alien_NT=res.NT.Alien_noise;
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;
Rx_LT=res.LT.Rx_signal;
Tx_LT=res.LT.Tx_signal;
tmp=find(Tx_LT>10.^-14);
NT_mask=10.^14.*ones(size(f));
NT_mask(tmp)=ones(size(tmp));

% Increase the noise with xtalk_margin 
noise_LT=Tot_LT.*10.^(xtalk_margin/10);
noise_NT=Tot_NT.*10.^(xtalk_margin/10);

SNR_LT=Rx_LT./noise_LT;
SNR_NT=Rx_NT./noise_NT;

H_LT=interp1(f,Rx_LT./Tx_NT,newf);
N_LT=interp1(f,noise_LT,newf);
H_NT=interp1(f,Rx_NT./Tx_LT,newf);
N_NT=interp1(f,noise_NT,newf);
PSD_NT=interp1(f,Tx_NT,newf);
PSD_LT=interp1(f,Tx_LT,newf);
Gap=10.^((refSNR+SNRloss-codingGain)/10);
E_total=10.^(Px./10)./df;
   
[b_LT, E_LT, g_LT, M_LT] =dmt(sqrt(abs(H_LT)),N_LT,E_total,Gap,PSD_NT,constellation);
[b_NT, E_NT, g_NT, M_NT] =dmt(sqrt(abs(H_NT)),N_NT,E_total,Gap,PSD_LT,constellation);
     
% And sum up the rates
index_LT=find((newf>=f_start_LT)&(newf<=f_stop_LT));
index_NT=find((newf>=f_start_NT)&(newf<=f_stop_NT));
Rate_LT=trapz(newf(index_LT),b_LT(index_LT))./1e6.*TimeDivision_up.*(1-efficiencyLoss);
Rate_NT=trapz(newf(index_NT),b_NT(index_NT))./1e6.*TimeDivision_down.*(1-efficiencyLoss);

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