function [dfseSNR,w_t,opt_delay]=dfsecolorsnr(l,p,nff,nbb,delay,Ex,noise);
% --------------------------------------------------------------
% [dfseSNR,w_t] = dfecolor(l,p,nff,nbb,delay,Ex,noise);
% l     = oversampling factor
% p     = pulse response, oversampled at l (size)
% nff   = number of feed-forward taps
% nbb   = number of feedback taps
% delay = delay of system <= nff+length of p - 2 - nbb
%         if delay = -1, then choose best delay
% Ex    = average energy of signals
% noise = noise autocorrelation vector (size l*nff)
% NOTE: noise is assumed to be stationary
%
% outputs:
% dfseSNR = equalizer SNR, unbiased in dB
% w_t = equalizer coefficients [w -b]
% opt_delay = optimal delay found if delay =-1 option used.
%             otherwise, returns delay value passed to function
% created 4/96;
% ---------------------------------------------------------------

size = length(p);
nu = ceil(size/l)-1;
p = [p zeros(1,(nu+1)*l-size)];

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

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

%form matrix P, vector channel matrix
P = zeros(nff*l+nbb,nff+nu);
for i=1:nff,
  P(((i-1)*l+1):(i*l),i:(i+nu)) = ptmp;
end

%precompute Rn matrix - constant for all delays
Rn  = zeros(nff*l+nbb);
Rn(1:nff*l,1:nff*l) = toeplitz(noise);

dfseSNR = -100;
P_init = P;
%loop over all possible delays
for d = delay,
    P = P_init;
    P(nff*l+1:nff*l+nbb,d+2:d+1+nbb) = eye(nbb);
    %P
    temp= zeros(1,nff+nu);
    temp(d+1)=1;
    %construct matrices
    Ry  = Ex*P*P' + Rn;
    Rxy = Ex*temp*P';
    new_w_t = Rxy*inv(Ry);
    sigma_dfse = Ex - real(new_w_t*Rxy');
    new_dfseSNR = 10*log10(Ex/sigma_dfse - 1);
    %save setting of this delay if best performance so far
    if new_dfseSNR >= dfseSNR
        w_t = new_w_t;
        dfseSNR = new_dfseSNR;
        opt_delay = d; 
    end
end

