function PSD=modelPSD_SDSL_sym(lcname,f)
%% ===========================================================================
% modelPSD_SDSL_sym - Calculates PSDs for symmetric SDSL coded 16 PAM 
%
% Parameter:    lcname    Linecode name
% 		f         Frequency
% Returns:      PSD       Resulting SDSL PSD
%
%
% Reference:    Draft ETSI TS 101524-2 V1.1.1 (2000-05)  
%% ===========================================================================

%% ===========================================================================
% Copyright (C) 2000 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                     All rights reserved. 
% Project       : FSAN duplex model
% Author(s)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%               : Gernot Schmid (schmid@ftw.at)
%
%% ===========================================================================
% Change History
%      2000-03-23 (GS) Created
%      2000-04-12 (GS) Extensions for optionally applying a receive filter
%      2000-09-27 (GS) Cleaning up for release
%% ===========================================================================
global ex

lc=getList(ex.lclist,lcname);

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

% Find correct scaling factor K for considered bitrate in template
k=lc.param.pam.k;
tmpk=k(1,:);
tmpind=find(tmpk>brate); 
K=k(2,tmpind(1)-1);
 
fsym =  (brate+overh)/bits_per_sym;          % Transmitted symbol rate
f3dB = fsym*lc.param.pam.filterratio;        % 3dB cutoff of transmit butterworth filter


% Receive filter 3dB frequencies (optionally used in 'calcResultSDSL.m')
lc.param.pam.recfilter3dBf.dn=f3dB; % (assumed as transmit filter) 
lc.param.pam.recfilter3dBf.up=f3dB; % (assumed as transmit filter) 
ex.lclist   = setList(ex.lclist,lcname,lc);  % Add it to linecode list
 

% Set up MaskOffset
MOa=1+0.4.*((f3dB.*ones(size(f))-f)./(f3dB.*ones(size(f))));
MOb=ones(size(f));
ind3dB=find(f>=f3dB);
MaskOffset(1:ind3dB(1))=MOa(1:ind3dB(1));
MaskOffset(ind3dB(1)+1:length(f))=MOb(ind3dB(1)+1:length(f));

MaskOffset=0; % Deactivates Maskoffset if nominal PSDs should be used


PSDa= 1e3.*((K.* 1./(fsym.*R).*(sinc(f./(fsym))).^2./(1 + (f./f3dB).^(2.*N))))...
      .*10.^(MaskOffset./10).*(f.^2./(f.^2+fc.^2)); % Main PSD in mW/Hz

PSDb=1e3.*ones(size(f)).*f.^(-1.5).*0.5683.*1e-4;   % Floor PSD in mW/Hz

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

%Resulting PSD:
PSD(1:intersect-1)=PSDa(1:intersect-1);
PSD(intersect:length(f))=PSDb(intersect:length(f));