function PSD=modelPSD_SDSL_sym_ver2(lcname,f)
%% ===========================================================================
% modelPSD_SDSL_sym - Calculates PSDs for symmetric SDSL coded 16 PAM 
%
% Parameter:    lcname    Linecode name
% 		f         Frequency
% Returns:      PSD       Resulting SDSL PSD
%
%
% References:
%   Draft ETSI TS 101524-2 V1.1.1 (2000-05)
%   
%   ETSI TR 101 830-2,  Spectral management on metallic access
%   networks; Part 2: Technical  methods and performance
%   evaluations, m01p20a5  
%% ===========================================================================

%% ===========================================================================
% Copyright (C) 2000-2009 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                     All rights reserved. 
% Project       : FTW's xDSLsimu
% Author(s)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%               : Gernot Schmid (schmid@ftw.at)
%
% CVS:       $Id: modelPSD_SDSL_sym_SpM.m,v 1.3 2005/01/04 10:32:28 tono Exp $
%% ===========================================================================
% Change History
%      2003-11-11 (ToNo) Created based on an earlier version
%      2004-01-14 (ToNo) Update according to ETSI SpM part 2 complete
%      2005-01-11 (ToNo) Renaming and integration into main code
%      2005-01-13 (ToNo) Set up to better follow the ETSI SpM part2
%      2008-08-28 (ToNo) Introduced consistent naming (pam -> mod)
%% ===========================================================================
global ex

lc=getList(ex.lclist,lcname);

bits_per_sym = lc.param.mod.bpsym;           % Payload bits per PAM symbol
brate        = lc.param.mod.brate.rate;      % Payload bit rate 
overh        = lc.param.mod.brate.ohead;     % Overhead
% fc           = lc.param.mod.highpass3dBf;    % Coupling transformer cutoff fq
R            = ex.param.Zterm; 	             % Line termination impedance
N            = lc.param.mod.filterorder;     % Order of transmit butterworth filter
floorf       = lc.param.mod.floorf;          % When the PSD start being flat


Kx           = 0.5683e-4;                    % [W/Hz]
fl           = 5e3;                          % [Hz]

% Find correct scaling factor K for considered bitrate in template
k=lc.param.mod.k;
tmpk=k(1,:);
tmpind=find(tmpk>brate);
K=k(2,tmpind(1)-1);

fhfactor     = 0.5;
fxfactor     = 2;

fsym         = (brate+overh)/bits_per_sym;   % Transmitted symbol rate
fx           = fsym*fxfactor;
fh           = fx*fhfactor;


P1 = 1e3* K/(fx*R) * sinc(f/fx).^2 .* (1 ./ (1+(fl./f).^2)) ...
     .* (1 ./ (1+(f./fh).^(2*N)));           % Main PSD in mW/Hz

P2 = 1e3 * Kx * f.^(-1.5);                   % Floor PSD in mW/Hz

% Find the best intersect point in frequency vector
indint=find((P2>=P1));
dind=diff(indint);
tmp1=find(dind~=1)+1;
intersect=indint(tmp1);

%Resulting PSD:
PSD(1:intersect-1)=P1(1:intersect-1);
PSD(intersect:length(f))=P2(intersect:length(f));
PSD(f>floorf)=10^-11;