function G=sqrcPSD(f,alpha,fs,fc,level,expCutoflevel)
%sqrcPSD - Squared Raised Cosine PSD mask  
%
%  G=sqrcPSD(f,alpha,fs,fc)
%
% Parameter:    f         frequency axis
% Parameter:    alpha     excess bandwidth (0<=alpha<=1)
% Parameter:    fs        vector of sampling freqeuncies
% Parameter:    fc        vector of carrier frequencies
% Parameter:    level     vector of power levels (at fc) (in db)
% Parameter:    expCutoflevel  exponent of the cut of level
%                               
% Example(s):
%   G=sqrcPSD((1:10000).*1e3,0.2,[2e6 4e6],[1.3e6 5e6]);
%
%   This function will generate a PSD mask that is square-root raised-cosine shaped 
%   for each index in the fs, fc vectors
%
% 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)
%               : 
%
% CVS:       $Id: sqrcPSD.m,v 1.5 2000/04/06 12:35:29 uvan Exp $
%% ===========================================================================
% Change History
%      1999-08-16 (DaB)  Created
%      2000-04-03 (UvAn) Possible to change the exponent of the cutoflevel with a parameter 
%% ===========================================================================

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

if nargin<5,
    level=ones(1,length(fc));
else
    level=10.^(level/10);
end

if nargin<6,
    expCutoflevel = -12;
end

G=zeros(1,length(f));
for a=1:length(fs),
   f_carr=fc(a);
   T=1./fs(a);
   tail_max=(1/(2*T).*(1+alpha))+f_carr;
   tail_min=(1/(2*T).*(1-alpha))+f_carr;
   f_tail=find((f>=tail_min)&(f<=tail_max));
   if alpha>0,
      tail=sqrt(0.5*(1+cos(pi*T/alpha.*(f(f_tail)-(1-alpha)./(2*T)-f_carr))));
   else
      tail=0;
   end;
   
   head_min=-(1/(2*T).*(1+alpha))+f_carr;
   head_max=-(1/(2*T).*(1-alpha))+f_carr;
   f_head=find((f>=head_min)&(f<=head_max));
   if alpha>0,
      head=sqrt(0.5*(1+cos(pi*T/alpha.*(f(f_head)+(1-alpha)./(2*T)-f_carr))));
   else
      head=0;
   end;
   
   f_main=find((f>head_max)&(f<tail_min)); 
   main=ones(1,length(f_main));

   G(f_head)=level(a).*head;
   G(f_main)=level(a).*main;
   G(f_tail)=level(a).*tail;
end;
cutoflevel=10^expCutoflevel;
G(find(G<=cutoflevel))=cutoflevel.*ones(size(find(G<=cutoflevel)));