function [f_opt,b_opt,delay,mmse]=dfe_tap_estimator(channel,Lf,Lb,lambda,P,bias);

% [f_opt,b_opt,delay,mmse]=dfe_tap_estimator(C,Lf,Lb,lambda,P,bias)
%
% This function performs the matrix based MMSE estimates of the DFE taps
% channel - Channel impulse response (should be a row vector)
% Lf - Number of FFF taps
% Lb - Number of FBF taps
% lambda - (Noise variance)/(Signal Variance)
% P - oversampling factor (either baud spaced P=1/ fractionally spaced P=2)
% bias  - 0=unbiased, 1=biased
% f_opt - optimal FFF filter taps
% b_opt - optimal FBF filter taps
% delay - optimal delay 
% mmse  - mimimum mean-squared error

[h]=channel;  % Initialize channel parameters;
Lh=length(h); % length of channel

cm = convmtx(h,Lf)'; % Find convolution matrix

if (P==2) %fractionally-spaced tap optimization

cm1=cm(1:2:Lf+Lh-1,:); % Down sampled convolution matrix, by T/2
[delta1,junk]=size(cm1);
cm2=cm(2:2:Lf+Lh-1,:); 
[delta2,junk]=size(cm2);

I_delta1=eye(delta1);
I_delta2=eye(delta2);

f1=zeros(Lf,delta1);
b1=zeros(Lb,delta1);
for delay1=1:delta1, % Delay optimization

I_delta=I_delta1;
I_delta(delay1+1:min(delta1,delay1+Lb),delay1+1:min(delta1,delay1+Lb))=zeros(min(Lb,delta1-delay1));
m1=inv(cm1'*I_delta*cm1+lambda*eye(Lf))*cm1'*I_delta ;
temp=I_delta*cm1*m1;
mse1(delay1)=1-temp(delay1,delay1); % MSE calculation
f1(:,delay1)=m1(:,delay1); % Feedforward equalizer taps
temp=(I_delta1-I_delta)*cm1*f1(:,delay1);
b1(1:min(Lb,delta1-delay1),delay1)=temp(delay1+1:min(delta1,delay1+Lb)); 

end;

for delay2=1:delta2,

I_delta=I_delta2;
I_delta(delay2+1:min(delta2,delay2+Lb),delay2+1:min(delta2,delay2+Lb))=zeros(min(Lb,delta2-delay2));
m2=inv(cm2'*I_delta*cm2+lambda*eye(Lf))*cm2'*I_delta ;
temp=I_delta*cm2*m2;
mse2(delay2)=1-temp(delay2,delay2);
f2(:,delay2)=m2(:,delay2);
temp=(I_delta2-I_delta)*cm2*f2(:,delay2);
b2(1:min(Lb,delta2-delay2),delay2)=temp(delay2+1:min(delta2,delay2+Lb));

end;

mse=zeros(length(mse1)+length(mse2),1);
mse(1:2:length(mse))=mse1;
mse(2:2:length(mse))=mse2;

[mmse,index]=min(mse);  % Delay optimization of mse

if (rem(index,2)~=0)
f_opt=f1(:,ceil(index/2)); % FFF taps
b_opt=b1(:,ceil(index/2)); % FBF taps
else
f_opt=f2(:,index/2);
b_opt=b2(:,index/2);
end;

delay=index;

elseif (P==1) %baud-spaced tap optimization

[delta1,junk]=size(cm);
I_delta1=eye(delta1);

f1=zeros(Lf,delta1);
b1=zeros(Lb,delta1);
for delay1=1:delta1, % Delay optimization

I_delta=I_delta1;
I_delta(delay1+1:min(delta1,delay1+Lb),delay1+1:min(delta1,delay1+Lb))=zeros(min(Lb,delta1-delay1));
m1=inv(cm'*I_delta*cm+lambda*eye(Lf))*cm'*I_delta ;
temp=I_delta*cm*m1;
mse1(delay1)=1-temp(delay1,delay1); % MSE calculation
f1(:,delay1)=m1(:,delay1); % Feedforward equalizer taps
temp=(I_delta1-I_delta)*cm*f1(:,delay1);
b1(1:min(Lb,delta1-delay1),delay1)=temp(delay1+1:min(delta1,delay1+Lb)); 

end;
[mmse,index]=min(mse1);  % Delay optimization of mse
f_opt=f1(:,(index)); % FFF taps
b_opt=b1(:,(index)); % FBF taps
delay=index;
end;

if (bias==0) % unbias the DFE taps
temp=conv(h,f_opt);
f_opt=f_opt/temp(delay);
b_opt=b_opt/temp(delay);
mmse=(mmse-(1-temp(delay))^2)/(temp(delay)^2);
end;









