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

%% ===========================================================================
% Copyright (C) 1998,99 by Telia Research AB, Lulea, Sweden; All rights reserved. 
% Project       : FSAN duplex model
% Author(s)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%               : Daniel Bengtsson (Daniel.J.Bengtsson@Telia.se)
%
% CVS:       $Id: calcPSD.m,v 1.9 2000/04/26 15:41:21 tono Exp $
%% ===========================================================================
% 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
%% ===========================================================================

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

if nargin<4,
    low=min(data(2:2:length(data)));
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'
                    tmp=1./((min(log10(f(x))))-(max(log10(f(x))))).* ...
                        (log10(f(x))-min(log10(f(x))));
                otherwise
                    fprintf('Undefined Interpolation');
                    break;
            end;
        end;
    end;  
    PSD(x)=level(index-1)+(level(index-1)-level(index))*tmp;   
end;

PSD=10.^(PSD/10);