function [SumSNR]=calcFoldedSNR(SNR,f,fsym,modulation);
%% ===========================================================================
% calc_fold_SNR - Calculates the folded SNR for DFE-based
%                 PAM/QAM-margin calculation (applicable even for 
%                 non equally spaced f-axis)
%
% Parameter:    SNR   	Signal to Noise ratio as simple Rx/Noise expression
% Parameter:    f	Frequency axis
% Parameter:    fsym	Symbol transmission rate
% Parameter:    modulation  'PAM', 'CAP', or 'QAM', where 'PAM' is default
%
% Returns:      SumSNR	Folded SNR in dB
%
% Reference: 	T1E1.4/2000-002R2 DRAFT ANSI T1.XXX-2000
%
%% ===========================================================================

%% ===========================================================================
% Copyright (C) 2000-2003 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                     All rights reserved. 
% Project       : FTW's xDSLsimu
% Author(s)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%               : Gernot Schmid (schmid@ftw.at)
%
%% ===========================================================================
% Change History
%      2000-05-30 (GS)   Created
%      2000-06-20 (ToNo) Cleaning up  
%      2000-07-11 (ToNo) Made it use the new way of interpolation
%      2003-11-26 (ToNo) Added QAM modulation
%      2003-12-19 (ToNo) Removed unneeded call to conj
%% ===========================================================================

if nargin<4,
    modulation = 'PAM';
end

% Generate a new f scale     
deltaf=100;
fnew=[deltaf:deltaf:fsym+1-deltaf];
fnewinv=fsym-fnew;

% Get the folds
SNR0  =   interpolate(f,SNR,fnew);           % fold 0
SNR1  =   interpolate(f,SNR,fnew+fsym);      % fold 1
switch modulation
 case 'PAM'
  SNR2  = interpolate(f,SNR,fnewinv);        % fold -1
  SNR3  = interpolate(f,SNR,fnewinv+fsym);   % fold -2
 case {'QAM', 'CAP'}
  SNR2  = interpolate(f,SNR,fnew+fsym*2);    % fold 2
  SNR3  = interpolate(f,SNR,fnew+fsym*3);    % fold 3
 otherwise
  error(['calcFoldedSNR do not implement ' mod ' at the moment']);
end

% Sum the folds
folded_SNR=SNR0+SNR1+SNR2+SNR3;
% testfoldSNR =  interpolate(fnew,folded_SNR,f);
% plot(f,10*log10(testfoldSNR))
% plot(f,10*log10(interpolate(fnew,folded_SNR,f)))
% axis([0,1.1*fsym,0,70])

SumSNR=trapz(fnew,10*log10(1+folded_SNR))/fsym;