function PSD=modelPSD_SDSL_asym(lcname,f,upordown)
%% ===========================================================================
%modelPSD_SDSL_asym - calculates PSDs for asymmetric SDSL coded 16 PAM 
%
% Parameter:    lcname    Linecode name
% 		f         Frequency
%               upordown  Switch determines if up- or downstream is considered 
% 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

fsym=(brate+overh)/bits_per_sym; 	     % Transmitted symbol rate

% Definition for upstream and downstream
switch upordown 
  case 1 % Downstream
     ktmpl   = lc.param.pam.k.dn;  % Scaling factor template (see lcDefsSDSL_asym)
     Ntmpl   = lc.param.pam.filterorder.dn; % Filter order template (see lcDefsSDSL_asym)
     fxtmpl  = [0 2048e3 2304e3 1e7; fsym 1370667 1541333 154133;];  
     f3dBtmpl= [0 2048e3 2304e3 1e7; lc.param.pam.filterratio*fsym 548267 578000 578000;];
     
     % Receive filter 3dB frequency (optionally used in 'calcResultSDSL.m'
     lc.param.pam.recfilter3dBf.dn=f3dBtmpl; % (assumed as transmit filter)
     
  case 2 % Upstream
     ktmpl   = lc.param.pam.k.up;  % scaling factor template (see lcDefsSDSL_asym)
     Ntmpl   = lc.param.pam.filterorder.up; % filter order template (see lcDefsSDSL_asym)
     fxtmpl  = [0 2048e3 2304e3 1e7; fsym 685333 770667 770667;]; 
     f3dBtmpl= [0 2048e3 2304e3 1e7; lc.param.pam.filterratio*fsym 342667 385333 385333;];
     
     %Receive filter 3dB frequency (optionally used in 'calcResultSDSL.m
     lc.param.pam.recfilter3dBf.up=f3dBtmpl; % (assumed as transmit filter) 

  otherwise
     error('Invalid third input argument in modelPSD_SDSL_asym');
  end
  
% Add receive filter characteristics to linecode list 
ex.lclist   = setList(ex.lclist,lcname,lc);

  
% Find correct values in templates (depending on bitrate)
tmpk=ktmpl(1,:);
tmpind=find(tmpk>brate);
K=ktmpl(2,tmpind(1)-1);      % Scaling factor

tmpN=Ntmpl(1,:);
tmpind=find(tmpN>brate);
N=Ntmpl(2,tmpind(1)-1);      % Transmit butterworth filter order

tmpfx=fxtmpl(1,:);
tmpind=find(tmpfx>brate);
fx=fxtmpl(2,tmpind(1)-1);    % Characteristic frequency to define shape of PSD 

tmpf3dB=f3dBtmpl(1,:);
tmpind=find(tmpf3dB>brate);
f3dB=f3dBtmpl(2,tmpind(1)-1); % 3dB cutoff frequency of transmit filter
 
% 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./(fx))).^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));