function Pout=modelPSD_HDSL2B1Q(f,pairs)
%% ===========================================================================
% modelPSD_HDSL2B1Q_v121 - PSD mask for HDSL 2B1Q
%
% Parameter:    f	Frequency axis
% Parameter:    pairs	No of pairs
% Returns:      Pout	PSD mask (in mW)
%
% Example(s):
%   
% Algorithmic details:
%   
% Reference:
%   Draft ETSI TR 101 830-2 V1.2.1 (2008-xx); pds:m06p06a02_SpM2_DR
%% ===========================================================================

%% ===========================================================================
% Copyright (C):
%       2008-2009 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                         All rights reserved.
% Project       : FTW's xDSLsimu
% Author(s)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%
% CVS:       $Id: $
%% ===========================================================================
% Change History
%      2008-03-17 (ToNo) Created based on "modelPSD_HDSL2B1Q"
%% ===========================================================================

if nargin<2,
    pairs = 2;
end

switch pairs
 case 1
  pvec = [1160e3, 3e3, 0.42*1160e3, 3, Inf, 1, 1.4662, 14, -133];
 case 2
  pvec = [584e3,  3e3, 0.68*584e3,  4, Inf, 1, 1.1965, 14, -133];
 case 21
  pvec = [584e3,  3e3, 0.50*584e3,  3, Inf, 1, 1.3501, 14, -133];
 case 22
  pvec = [584e3,  3e3, 0.68*584e3,  4, 1.5*584e3, 2, 1.1965, 14, -133  ];
 case 3
  pvec = [392e3,  3e3, 0.50*392e3,  3, Inf, 1, 1.3642, 14, -117.5];
 otherwise
  error('The number of pairs in modelPSD_HDSL2B1Q should be in the range 1-3')
end

fx     = pvec(1);
fl     = pvec(2);
fh1    = pvec(3);
Nh1    = pvec(4);
fh2    = pvec(5);
Nh2    = pvec(6);
qN     = pvec(7);
Ptot   = 10^(pvec(8)/10);
Pfloor = 10^(pvec(9)/10);


P1 = Ptot * ( 2 * qN /fx ) * sinc(f/fx).^2 .* (1 ./ (1+(fl./f).^2)) ...
     .* (1 ./ (1+(f./fh1).^(2*Nh1))) .* (1 ./ (1+(f./fh2).^(2*Nh2)));

Pout = max(P1,Pfloor);