function [F_opt, B_opt, Delay, Mse] = mimodfe_tap_est(C, Lf, Lb, M, P, Gamma, method, bias, samp);

% [F_opt, B_opt, Delay, Mse] = mimodfe_tap_est(C, Lf, Lb, M, P, Gamma, method, bias, samp);
%
% This function computes the MMSE-DFE coefficients for a MIMO channel with
% M transmitter antennas and P receiver antennas, P >= M by optimizing over
% all possible detection delays.
%
% Inputs: C - Chan. imp. response (MLc x P matrix), {\bf r}={\bf C}^{T}{\bf s}
%              _                                    _         _          _
%             | c11(0)    c12(0)    ...    c1P(0)    |       | s1(n)      |
%             |   :         :        :       :       |       |   :        |
%             | c11(Lc-1) c12(Lc-1) ...    c1P(Lc-1) |       | s1(n-Lc+1) |  
%         C = |   :         :        :       :       | , s = |   :        |
%             | cM1(0)    cM2(0)    ...    cMP(0)    |       | sM(n)      |
%             |   :         :        :       :       |       |   :        |
%             | cM1(Lc-1) cM2(Lc-1) ...    cMP(Lc-1) |       | sM(n-Lc+1) |
%             |_                                    _|       |_          _|
%
%         Lf - # of feedforward filter taps/receiver antenna
%         Lb - # of strictly causal feedback filter taps/transmitter antenna
%         M  - # transmit antennas
%         P  - # receive antennas
%         Gamma - the noise correlation matrix normalized by source variance
%                 od size (PLf x PLf)
%         method - 'joint'     : FBF cancels post-cursor ISI+MAI
%                - 'joint_sic' : FBF cancels post-cursor ISI+MAI and cursor
%                                  MAI of detected users.
%                - 'joint_sicrev': reversed order of detection
%         bias   - 0 = unbiased, 1 = biased.
%         samp   - oversampling factor (always use 1).
%
% Outputs: F_opt - Feedforward filter coefficients (PLf x M matrix), y=F^{T}r
%          B_opt - Feedforward filter coefficients (M(Lb+1) x M matrix),
%          Delay - vector of optimal delays, (M x 1)
%          Mse   -  vectore of mean squared erros, (M x 1).
%

% Lh = Lf+Lc-1
% Phi - channel convolution matrix (MLh x PLf)

Lc=size(C,1)/M;
Lh=Lc+Lf-1;
Phi=zeros(M*Lh,P*Lf); 

% Create channel convolution matrix

for m=1:M
for p=1:P
     tmp=convmtx(C((m-1)*Lc+1:m*Lc,p),Lf);
     Phi((m-1)*Lh+1:m*Lh,(p-1)*Lf+1:p*Lf)=tmp;
end;
end;

% Determine DFE tap coefficients

for delay=1:Lh-Lb, %delta+1
%delay=1; %fixed delay delta+1

b_opt=zeros(M*(Lb+1),M);
e=zeros(M*Lh,1); % to be a unit vector

V=eye(M*Lh,M*Lh); % to be a (MLh x MLh) matrix

for m=1:M,
V((m-1)*Lh+delay+1:(m-1)*Lh+delay+Lb,(m-1)*Lh+delay+1:(m-1)*Lh+delay+Lb)=zeros(Lb,Lb); 
end;

%%%%
if strcmp(method,'joint')

for m=1:M

Vm=V;

e_m=e;
e_m((m-1)*Lh+delay)=1;

f = pinv((Vm*Phi)'*(Vm*Phi)+Gamma)*(Vm*Phi)'*e_m;

f_opt(:,m)=f;%(:,delay);
b_tmp=Phi*f_opt(:,m);

for k=1:M,
b_opt(k*(Lb+1)-Lb:k*(Lb+1),m)=b_tmp((k-1)*Lh+delay:(k-1)*Lh+delay+Lb); 
end;

if (bias==0)
     f_opt(:,m)=f_opt(:,m)/b_opt(m*(Lb+1)-Lb,m);
     b_opt(:,m)=b_opt(:,m)/b_opt(m*(Lb+1)-Lb,m);
end;
mse(m,delay)=norm(Vm*Phi*f_opt(:,m)-e_m).^(2)+f_opt(:,m)'*Gamma*f_opt(:,m);
end;

end;

%%%%
if strcmp(method,'joint_sic')

Vm=V;

for m=1:M

if (m>1)
Vm((m-2)*Lh+delay,(m-2)*Lh+delay)=0;
end

e_m=e;
e_m((m-1)*Lh+delay)=1;

f = pinv((Vm*Phi)'*(Vm*Phi)+Gamma)*(Vm*Phi)'*e_m;

f_opt(:,m)=f;%(:,delay);
b_tmp=Phi*f_opt(:,m);

for k=1:M,
b_opt(k*(Lb+1)-Lb:k*(Lb+1),m)=b_tmp((k-1)*Lh+delay:(k-1)*Lh+delay+Lb); 
end;

if (bias==0)
     f_opt(:,m)=f_opt(:,m)/b_opt(m*(Lb+1)-Lb,m);
     b_opt(:,m)=b_opt(:,m)/b_opt(m*(Lb+1)-Lb,m);
end;
mse(m,delay)=norm(Vm*Phi*f_opt(:,m)-e_m).^(2)+f_opt(:,m)'*Gamma*f_opt(:,m);
end;

end;

%%%%
if strcmp(method,'joint_sicrev')

Vm=V;

for m=M:-1:1

if (m<M)
Vm(m*Lh+delay,m*Lh+delay)=0;
end

e_m=e;
e_m((m-1)*Lh+delay)=1;

f = pinv((Vm*Phi)'*(Vm*Phi)+Gamma)*(Vm*Phi)'*e_m;

f_opt(:,m)=f;%(:,delay);
b_tmp=Phi*f_opt(:,m);

for k=1:M,
b_opt(k*(Lb+1)-Lb:k*(Lb+1),m)=b_tmp((k-1)*Lh+delay:(k-1)*Lh+delay+Lb); 
end;

if (bias==0)
     f_opt(:,m)=f_opt(:,m)/b_opt(m*(Lb+1)-Lb,m);
     b_opt(:,m)=b_opt(:,m)/b_opt(m*(Lb+1)-Lb,m);
end;
mse(m,delay)=norm(Vm*Phi*f_opt(:,m)-e_m).^(2)+f_opt(:,m)'*Gamma*f_opt(:,m);
end;

end;

F(:,:,delay)=f_opt;
B(:,:,delay)=b_opt;

end

asnr=-10*log10(sum(mse,1));

[val,Delay]=max(asnr);

Mse=mse(:,Delay);
F_opt=F(:,:,Delay);
B_opt=B(:,:,Delay);












