% function [dfseSNR,W,b, optdelay]=dfeRAKE(l,h,nff,nbb,delay,Ex,noise);
% DFE design program for RAKE receiver
%
%
% Inputs
% l     = oversampling factor
% L is derived as No. of fingers in RAKE (number of rows in h)
% h     = pulse response matrix, oversampled at l (size), 
%         each row corresponding to a diversity path
% nff   = number of feedforward taps for each RAKE finger
% nbb   = number of feedback taps 
% delay = delay of system <= nff+length of p - 2 - nbb
% Ex    = average energy of signals
% noise = noise autocorrelation vector (size L x l*nff)
% NOTE: noise is assumed to be stationary, but may be spatially correlated
%
% outputs:
% dfseSNR = equalizer SNR, unbiased in dB
% W is feedforward section
% b is feedback filter
% optdelay is optimum delay
% ------------------------------
function [dfseSNR,W,b, opt_delay]=dfeRAKE(l,h,nff,nbb,delay,Ex,noise);

siz = size(h,2);
L=size(h,1);
nu = ceil(siz/l)-1;
h = [h zeros(L,(nu+1)*l-siz)];

% error checks
if nff<=0
  error('number of feedforward taps must be > 0');
end
if nbb<0
  error('number of feedback taps must be >= 0');
end
if delay > (nff+nu-1)
  error('delay must be <=  (nff+(length of p)-2)');
elseif delay < -1 
  error('delay must be >= -1');
elseif delay == -1
    delay = 0:1:nff+nu-1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%if length(noise) ~= L*l*nff
%   error('Length of noise autocorrelation vector must be L x(l*nff');
%end
if size(noise,2) ~= l*nff | size(noise,1) ~= L
   error('Size of noise autocorrelation matrix must be L x l*nff ');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%form ptmp = [p_0 p_1 ... p_nu] where p_i=[p(i*l) p(i*l-1)... p((i-1)*l+1)
for m=1:L
   ptmp((m-1)*l+1:m*l,1) = [h(m,1); zeros((l-1),1)];
end
for k=1:nu
   for m=1:L
      ptmp((m-1)*l+1:m*l,k+1) = conj((h(m,k*l+1:-1:(k-1)*l+2))');
   end
end

dfseSNR = -100;
%loop over all possible delays
for d = delay
    %% COMPUTE EXACT number of feedback taps WE NEED to 
    %% consider given delay d
    nbb_used = min(nff+nu-1-d,nbb);
    %form matrix P, vector channel matrix
    % form matrix P, vector channel matrix
    P = zeros(nff*l*L+nbb_used,nff+nu);

    %First construct the P matrix as in MMSE-LE
    for k=1:nff,
      P(((k-1)*l*L+1):(k*l*L),k:(k+nu)) = ptmp;
    end

    %A dd in part needed for the feedback
    P(nff*l*L+1:nff*l*L+nbb_used,d+2:d+1+nbb_used) = eye(nbb);
    temp= zeros(1,nff+nu);
    temp(d+1)=1;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Rn  = zeros(nff*l*L+nbb_used);

    for i = 1:L
       n_t = toeplitz(noise(i,:));
    
        for j = 1:l:l*nff
            for k = 1:l:l*nff
                Rn((i-1)*l+(j-1)*L+1:i*l+(j-1)*L, (i-1)*l+(k-1)*L+1:i*l+(k-1)*L) = n_t(j:j+l-1,k:k+l-1);
              end
        end
    end
    Rn;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    Ry  = Ex*P*P' + Rn;
    Rxy = Ex*temp*P';
%
    new_w_t = Rxy*inv(Ry);
    %Reshape the w_t matrix into the RAKE filter bank and feedback matrices
    ww=reshape(new_w_t(1:nff*l*L),l*L,nff);
    for m=1:L
      Wnew(m,:)=reshape(ww((m-1)*l+1:m*l,1:nff),1,l*nff);
    end
    bnew=-new_w_t(nff*l*L+1:nff*l*L+nbb_used);
    sigma_dfse = Ex - new_w_t*Rxy';
    new_dfseSNR = 10*log10(Ex/sigma_dfse - 1);

% save setting of this delay if best performance yet
    if new_dfseSNR >= dfseSNR
        W=Wnew;
        b=bnew;
        dfseSNR = new_dfseSNR;
        opt_delay = d; 
        nbb_final = nbb_used;
    end
end

if nbb_final < nbb
    warning(strcat('For OPTIMAL FB filter settings N_bb=',num2str(nbb_final),' was used instead of the N_bb=',num2str(nbb), ' inputed'));
end
