function [mse,rho,alpha]=trdfe_mse(C,f1,b1,delay1,f2,b2,delay2,lambda);

% [mse,rho,alpha]=trdfe_mse(C,f1,b1,delay1,f2,b2,delay2,lambda);
%
% This function calculates the mse for the LC-BiDFE consisting of
% normal-mode and time-reversal mode unbiased DFEs. Channel and DFE are 
% assumed to be symbol-spaced. The optimal weighting factor alpha is  
% calculated
%
% C - Channel impulse response (should be a row vector)
% f1 - optimal FFF filter taps, normal mode - column vectors
% b1 - optimal FBF filter taps, normal mode - column vectors
% delay1 - system delay, normal mode
% f2 - optimal FFF filter taps for time-reversal mode
% b3 - optimal FBF filter taps for time-reversal mode
% delay2 - time-reversal mode system delay 
% lambda - noise variance/signal variance
%
% mse - overall mean-squared error of the LC-BiDFE
% rho - the cross-correlation coefficient for the noise sequences
% alpha - the optimal weighting factor for the linear combining scheme

Lc=length(C); % length of channel
Lf1=length(f1);
Lb1=length(b1);
Lf2=length(f2);
Lb2=length(b2);

h1=conv(C,transpose(f1));
Lh1=length(h1);
h2=conv(fliplr(C),transpose(f2));
Lh2=length(h2);

% Calculate mse1, mse2 and rho
h1(delay1:delay1+Lb1)=zeros(1,Lb1+1);
h1=h1(1:Lh1);
h2(delay2:delay2+Lb2)=zeros(1,Lb2+1);
h2=h2(1:Lh2);

mse1=h1*h1'+f1'*f1*lambda;
mse2=h2*h2'+f2'*f2*lambda;

f2_hat=zeros(1,Lh1+Lh2+Lf1-2);
h2_hat=zeros(1,Lh1+Lh2+Lh1-2);
f2_hat(Lh1:Lh1+Lf2-1)=f2';
h2_hat(Lh1:Lh1+Lh2-1)=h2';

xcorr=f2_hat(Lf1+delay1+delay2-2:-1:delay1+delay2-1)*f1*lambda;
xcorr=xcorr+h2_hat(Lh1+delay1+delay2-2:-1:delay1+delay2-1)*h1';
rho=xcorr/sqrt(mse1*mse2);

alpha=(mse2-rho*sqrt(mse1*mse2))/(mse1+mse2-2*rho*sqrt(mse1*mse2));

mse=alpha^2*mse1+(1-alpha)^2*mse2+2*alpha*(1-alpha)*rho*sqrt(mse1*mse2);




