function PSD=calcPSD(data,method,f,low,fbreak)
%% ===========================================================================
%calcPSD - From a PSD description generate a PSD mask
%
% Parameter:    data	PSD mask description
% Parameter:    method	Interpolation method: 'Linear','Log-Linear', 'VDSL2')
% Parameter:    f       Frequency axis to use
% Parameter:    low     Value for no signal
% Parameter:    fbreak  Break frequency for VDSL2 interpolation change
% Returns:      PSD	The resulting PSD mask
%
% Example(s):
%   PSD=calcPSD(tfplan.upstream,tfplan.interpolation,ex.param.frequency.f);
%
% Algorithmic details:
%
% Reference:
%% ===========================================================================

%% ===========================================================================
% Copyright (C):                                        
%       1998-1999 by Telia Research AB, Lulea, Sweden;                
%       2000-2009 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                         All rights reserved.
% Project       : FTW's xDSLsimu
% Author(s)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%               : Daniel Bengtsson (Daniel.J.Bengtsson@Telia.se)
%
% CVS:       $Id: calcPSD.m 752 2009-01-02 13:03:52Z tono $
%% ===========================================================================
% Change History
%      1998-11-09 (DaB)  Created
%      1998-11-11 (ToNo) Added header and comments
%      1998-01-15 (DaB/ToNo) Fixed the glitch at last up/down switch
%      1998-01-17 (ToNo) Added new parameters for special f axis and low values
%      1999-09-30 (DaB)  Rewritten for new ex struct
%      2000-04-03 (UvAn) Removed global ex; frequency axis passed in function
%      2002-07-30 (ToNo) Removed a duplication of log10 calculations = speed!
%      2008-10-01 (ToNo) Added VDSL2 type of mixed interpolation
%% ===========================================================================

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

if (nargin<4) || (low == 0)
    low=min(data(2:2:length(data)));
end
   
if nargin<5,
    fbreak = 3750e3-175e3; % VDSL2 Annex A
end
   
data=[f(1) low data f(length(f)) low]; % Fix glitch in last band
freq=data(1:2:length(data));
level=data(2:2:length(data));

PSD=data(length(data)) * ones(size(f));   	% High frequencies
tmp=find(f<=freq(1));                   
PSD(tmp)=level(1).*ones(size(tmp));     	% Low frequencies

for index=2:length(freq),               	% Set up data 
    x=find((f>=freq(index-1))&(f<=freq(index)));
    if isempty(x) || (length(x)==1),
        tmp=0;
    else   
        if (min(f(x))==max(f(x))),
            tmp=0;   
        else  	
            switch method
              case 'Linear'
                tmp=1./((min(f(x)))-max(f(x))).*(f(x)-min(f(x)));
              case 'Log-Linear'
                log10f=log10(f(x));
                tmp=1./((min(log10f))-(max(log10f))).* (log10f-min(log10f));
              case 'VDSL2'
                % Should actually check if fbreak is somewhere in
                % the range, but this makes it just too complicated
                if fbreak < freq(index)
                    log10f=log10(f(x));
                    tmp=1./((min(log10f))-(max(log10f))).* (log10f-min(log10f));
                else
                    tmp=1./((min(f(x)))-max(f(x))).*(f(x)-min(f(x)));
                end;
              otherwise
                fprintf('Undefined Interpolation');
                break;
            end;
        end;
    end;  
    PSD(x)=level(index-1)+(level(index-1)-level(index))*tmp;   
end;

PSD=10.^(PSD/10);