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:
%    FTW's xDSLsimu manual
%% ===========================================================================

%% ===========================================================================
% Copyright (C):                                        
%       1998-1999 by Telia Research AB, Lulea, Sweden;                
%       2000-2003 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,v 3.2 2003/03/19 11:14:41 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
%      2002-07-30 (ToNo) Removed a duplication of log10 calculations = speed!
%% ===========================================================================

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'
                    log10f=log10(f(x));
                    tmp=1./((min(log10f))-(max(log10f))).* ...
                        (log10f-min(log10f));
                otherwise
                    fprintf('Undefined Interpolation');
                    break;
            end;
        end;
    end;  
    PSD(x)=level(index-1)+(level(index-1)-level(index))*tmp;   
end;

PSD=10.^(PSD/10);