function [mse,alpha]=bidfe_mse_mimo(C1,F1,B1,delay1,method1,C2,F2,B2,delay2,method2, lambda,M,P);

% [mse,alpha]=bidfe_mse_mimo(C1,F1,B1,delay1,method1,C2,F2,B2,delay2,method2,lambda,M,P);
%
% This function calculates the mse for a linear combining of normal-mode
% and reverse-mode operation of an unbiased MIMO-DFE. Channel is assumed to be
% symbol-spaced, DFE is symbol-spaced. The optimal weighting factor alpha is 
% calculated
%
% C1 - Channel impulse response matrix (see dfe_tap_estmimo.m for details)
% F1 - optimal FFF filter taps, normal mode - matrix (PLf x M)
% B1 - optimal FBF filter taps, normal mode - matrix (M(Lb+1) x M)
% delay1 - system delay, normal mode
% method1 - type of MIMO DFE ('joint', 'joint_sic', 'joint_sicrev')
% C2 - Time-reversed channel response (MLc x P)
% F2 - optimal FFF filter taps, tr mode - matrix (PLf x M)
% B2 - optimal FBF filter taps, tr mode - matrix (M(Lb+1) x M)
% delay2 - tr mode system delay 
% method2 - type of MIMO trDFE ('joint', 'joint_sic', 'joint_sicrev')
% lambda - noise variance/signal variance
% M - # tx antennas
% P - # rx antennas

% Outputs: mse - a vector of mse values
%          alpha - a vector of weighting coefficients

Lc=size(C1,1)/M; % length of channel
Lf1=size(F1,1)/P;
Lb1=size(B1,1)/M-1;
Lf2=size(F2,1)/P;
Lb2=size(B2,1)/M-1;
Lh1=Lc+Lf1-1;
Lh2=Lc+Lf2-1;

% Generate channel convolution matrix

for m=1:M
for p=1:P
     tmp=convmtx(C1((m-1)*Lc+1:m*Lc,p),Lf1);
     Phi1((m-1)*Lh1+1:m*Lh1,(p-1)*Lf1+1:p*Lf1)=tmp;
     tmp=convmtx(C2((m-1)*Lc+1:m*Lc,p),Lf2);
     Phi2((m-1)*Lh2+1:m*Lh2,(p-1)*Lf2+1:p*Lf2)=tmp;
end;
end;

H1=Phi1*F1;
H2=Phi2*F2;
V1=eye(M*Lh1,M*Lh1); 
V2=eye(M*Lh2,M*Lh2);

% Generate V matrices 
for m=1:M,
V1((m-1)*Lh1+delay1+1:(m-1)*Lh1+delay1+Lb1,(m-1)*Lh1+delay1+1:(m-1)*Lh1+delay1+Lb1)=zeros(Lb1,Lb1);
V2((m-1)*Lh2+delay2+1:(m-1)*Lh2+delay2+Lb2,(m-1)*Lh2+delay2+1:(m-1)*Lh2+delay2+Lb2)=zeros(Lb2,Lb2);
end;

% Calculate mse1, mse2 and rho

G1=V1*H1;
G2=V2*H2;

for m=1:M

G1((m-1)*Lh1+delay1,m)=0;
G2((m-1)*Lh2+delay2,m)=0;

if strcmp(method1,'joint_sic')
G1(delay1:Lh1:(m-2)*Lh1+delay1,m)=0;
elseif strcmp(method1,'joint_sicrev')
G1((m-1)*Lh1+delay1:Lh1:M*Lh1,m)=0;
end;

if strcmp(method2,'joint_sic')
G2(delay2:Lh2:(m-2)*Lh2+delay2,m)=0;
elseif strcmp(method2,'joint_sicrev')
G2((m-1)*Lh2+delay2:Lh2:M*Lh2,m)=0;
end;

mse1(m)=norm(G1(:,m))^2+lambda*norm(F1(:,m))^2;
mse2(m)=norm(G2(:,m))^2+lambda*norm(F2(:,m))^2;

xcorr=0;
for k=1:M,
g2=zeros(1,Lh1+Lh2+Lh1-2);
g2(Lh1:Lh1+Lh2-1)=G2((k-1)*Lh2+1:k*Lh2,m)';
xcorr=xcorr+g2(Lh1+delay1+delay2-2:-1:delay1+delay2-1)*G1((k-1)*Lh1+1:k*Lh1,m);
end;
for k=1:P
f2=zeros(1,Lh1+Lh2+Lf1-2);
f2(Lh1:Lh1+Lf2-1)=F2((k-1)*Lf2+1:k*Lf2,m)';
xcorr=xcorr+f2(Lf1+delay1+delay2-2:-1:delay1+delay2-1)*F1((k-1)*Lf1+1:k*Lf1,m)*lambda;
end;

% Compute alpha, mse terms

alpha(m)=(mse2(m)-xcorr)/(mse1(m)+mse2(m)-2*xcorr);
mse(m)=alpha(m)^2*mse1(m)+(1-alpha(m))^2*mse2(m)+2*alpha(m)*(1-alpha(m))*xcorr;

end;






