function [ser,ser_tr,ser_soft,ser_soft_tr,ser_ideal,ser_ideal_tr,ser_bidfe,ser_soft_bidfe]=bidfe_ser_mimo(C1,F1,B1,delay1,method1,Mse1,C2,F2,B2,delay2,method2,Mse2,Alpha,snr,M,P,block_num,N);

% MIMO BiDFE, SER simulations; channel can be asymmetric
% [mse,alpha]=bidfe_mse_mimo(C1,F1,B1,delay1,method1,C2,F2,B2,delay2,method2,snr,M,P,block_num,N);
%
% This function calculates the mse for a diversity combining of normal-mode
% and reverse-mode operation of an unbiased MIMO-DFE. Channel is assumed to be
% asymmetric, 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')
% Mse1 - vector of MSE values at the output of the DFE
% 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')
% Mse2 - vector of MSE values at the output of the DFE
% Alpha - combining factor for BiDFE
% snr - SNR in dB
% M - # tx antennas
% P - # rx antennas
% block_num - number of test blocks
% N - length of data packet
% As of now 'joint' has not been implemented

% Outputs: ser - a vector of ser values
%          ser_tr - a vector of ser values of Time-reversal MIMO DFE
%          ser_soft - a vector of ser values of Soft MIMO DFE
%          ser_soft_tr - a vector of ser values of Soft Time-reversal MIMO DFE
%          ser_ideal - a vector of ser values of Ideal MIMO DFE
%          ser_ideal_tr - a vector of ser values of Ideal Time-reversal MIMO DFE
%          ser_bidfe - a vector of ser values of MIMO BiDFE
%          ser_soft_bidfe - a vector of ser values of Soft MIMO BiDFE

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;

tic

for iter=1:block_num

% terminal output
      fprintf('%4i ',iter);
      if rem(iter,20)==0
         fprintf('\n');
      end

% Generate Source Sequence
for m=1:M
symb(m,:)=pam(N,2,1);
end;

% Transmit through channel
rec_seq=zeros(P,N+Lc-1);
for p=1:P
for m=1:M
rec_seq(p,:)=rec_seq(p,:)+conv(symb(m,:),C1((m-1)*Lc+1:m*Lc,p));
end;
rec_seq(p,:)=rec_seq(p,:)+noise(length(rec_seq(p,:)),snr,1);
end;

% Feedforward filtering
filt_seq=zeros(M,N+Lc+Lf1-2);
for m=1:M
for p=1:P
filt_seq(m,:)=filt_seq(m,:)+conv(rec_seq(p,:),F1((p-1)*Lf1+1:p*Lf1,m));
end;
end;

filt_seq=filt_seq(:,delay1:end);

% Feedback filtering

% Ideal feedback
symb_hat_ideal(:,1:Lb1)=symb(:,1:Lb1);
for k=Lb1+1:N
if strcmp(method1,'joint_sic')
for m=1:M
y(m,k)=filt_seq(m,k);
for m1=1:M
y(m,k)=y(m,k)-symb(m1,k-1:-1:k-Lb1)*B1((m1-1)*(Lb1+1)+2:m1*(Lb1+1),m);
end;
if (m>1)
y(m,k)=y(m,k)-symb(1:m-1,k)'*B1(1:Lb1+1:(m-1)*(Lb1+1),m);
end;
symb_hat_ideal(m,k)=pamdet(y(m,k),2,1);
end;
elseif strcmp(method1,'joint_sicrev')
for m=M:-1:1
y(m,k)=filt_seq(m,k);
for m1=1:M
y(m,k)=y(m,k)-symb(m1,k-1:-1:k-Lb1)*B1((m1-1)*(Lb1+1)+2:m1*(Lb1+1),m);
end;
if (m<M)
y(m,k)=y(m,k)-symb(m+1:M,k)'*B1(m*(Lb1+1)+1:(Lb1+1):M*(Lb1+1),m);
end;
symb_hat_ideal(m,k)=pamdet(y(m,k),2,1);
end;
end;
end;

% Hard Decision feedback
symb_hat(:,1:Lb1)=symb(:,1:Lb1);
for k=Lb1+1:N
if strcmp(method1,'joint_sic')
for m=1:M
y(m,k)=filt_seq(m,k);
for m1=1:M
y(m,k)=y(m,k)-symb_hat(m1,k-1:-1:k-Lb1)*B1((m1-1)*(Lb1+1)+2:m1*(Lb1+1),m);
end;
if (m>1)
y(m,k)=y(m,k)-symb_hat(1:m-1,k)'*B1(1:Lb1+1:(m-1)*(Lb1+1),m);
end;
symb_hat(m,k)=pamdet(y(m,k),2,1);
end;
elseif strcmp(method1,'joint_sicrev')
for m=M:-1:1
y(m,k)=filt_seq(m,k);
for m1=1:M
y(m,k)=y(m,k)-symb_hat(m1,k-1:-1:k-Lb1)*B1((m1-1)*(Lb1+1)+2:m1*(Lb1+1),m);
end;
if (m<M)
y(m,k)=y(m,k)-symb_hat(m+1:M,k)'*B1(m*(Lb1+1)+1:(Lb1+1):M*(Lb1+1),m);
end;
symb_hat(m,k)=pamdet(y(m,k),2,1);
end;
end;
end;

% Soft Decision feedback

symb_hat_soft(:,1:Lb1)=symb(:,1:Lb1);
for k=Lb1+1:N
if strcmp(method1,'joint_sic')
for m=1:M
y_soft(m,k)=filt_seq(m,k);
for m1=1:M
y_soft(m,k)=y_soft(m,k)-symb_hat_soft(m1,k-1:-1:k-Lb1)*B1((m1-1)*(Lb1+1)+2:m1*(Lb1+1),m);
end;
if (m>1)
y_soft(m,k)=y_soft(m,k)-symb_hat_soft(1:m-1,k)'*B1(1:Lb1+1:(m-1)*(Lb1+1),m);
end;
symb_hat_soft(m,k)=tanh(y_soft(m,k)/sqrt(Mse1(m)));
end;
elseif strcmp(method1,'joint_sicrev')
for m=M:-1:1
y_soft(m,k)=filt_seq(m,k);
for m1=1:M
y_soft(m,k)=y_soft(m,k)-symb_hat_soft(m1,k-1:-1:k-Lb1)*B1((m1-1)*(Lb1+1)+2:m1*(Lb1+1),m);
end;
if (m<M)
y_soft(m,k)=y_soft(m,k)-symb_hat_soft(m+1:M,k)'*B1(m*(Lb1+1)+1:(Lb1+1):M*(Lb1+1),m);
end;
symb_hat_soft(m,k)=tanh(y_soft(m,k)/sqrt(Mse1(m)));
end;
end;
end;

%%%%%%% Time reversed decision feedback

%Time reversal
rec_seq_tr=fliplr(rec_seq);
symb_tr=fliplr(symb);

% Feedforward filtering
filt_seq_tr=zeros(M,N+Lc+Lf2-2);
for m=1:M
for p=1:P
filt_seq_tr(m,:)=filt_seq_tr(m,:)+conv(rec_seq_tr(p,:),F2((p-1)*Lf2+1:p*Lf2,m));
end;
end;

filt_seq_tr=filt_seq_tr(:,delay2:end);

% Feedback filtering

% Ideal feedback
symb_hat_ideal_tr(:,1:Lb2)=symb_tr(:,1:Lb2);
for k=Lb2+1:N
if strcmp(method2,'joint_sic')
for m=1:M
y_tr(m,k)=filt_seq_tr(m,k);
for m1=1:M
y_tr(m,k)=y_tr(m,k)-symb_tr(m1,k-1:-1:k-Lb2)*B2((m1-1)*(Lb2+1)+2:m1*(Lb2+1),m);
end;
if (m>1)
y_tr(m,k)=y_tr(m,k)-symb_tr(1:m-1,k)'*B2(1:Lb2+1:(m-1)*(Lb2+1),m);
end;
symb_hat_ideal_tr(m,k)=pamdet(y_tr(m,k),2,1);
end;
elseif strcmp(method2,'joint_sicrev')
for m=M:-1:1
y_tr(m,k)=filt_seq_tr(m,k);
for m1=1:M
y_tr(m,k)=y_tr(m,k)-symb_tr(m1,k-1:-1:k-Lb2)*B2((m1-1)*(Lb2+1)+2:m1*(Lb2+1),m);
end;
if (m<M)
y_tr(m,k)=y_tr(m,k)-symb_tr(m+1:M,k)'*B2(m*(Lb2+1)+1:(Lb2+1):end,m);
end;
symb_hat_ideal_tr(m,k)=pamdet(y_tr(m,k),2,1);
end;
end;
end;

% Hard Decision feedback

symb_hat_tr(:,1:Lb2)=symb_tr(:,1:Lb2);
for k=Lb2+1:N
if strcmp(method2,'joint_sic')
for m=1:M
y_tr(m,k)=filt_seq_tr(m,k);
for m1=1:M
y_tr(m,k)=y_tr(m,k)-symb_hat_tr(m1,k-1:-1:k-Lb2)*B2((m1-1)*(Lb2+1)+2:m1*(Lb2+1),m);
end;
if (m>1)
y_tr(m,k)=y_tr(m,k)-symb_hat_tr(1:m-1,k)'*B2(1:Lb2+1:(m-1)*(Lb2+1),m);
end;
symb_hat_tr(m,k)=pamdet(y_tr(m,k),2,1);
end;
elseif strcmp(method2,'joint_sicrev')
for m=M:-1:1
y_tr(m,k)=filt_seq_tr(m,k);
for m1=1:M
y_tr(m,k)=y_tr(m,k)-symb_hat_tr(m1,k-1:-1:k-Lb2)*B2((m1-1)*(Lb2+1)+2:m1*(Lb2+1),m);
end;
if (m<M)
y_tr(m,k)=y_tr(m,k)-symb_hat_tr(m+1:M,k)'*B2(m*(Lb2+1)+1:(Lb2+1):end,m);
end;
symb_hat_tr(m,k)=pamdet(y_tr(m,k),2,1);
end;
end;
end;

% Soft Decision feedback
symb_hat_soft_tr(:,1:Lb2)=symb_tr(:,1:Lb2);
for k=Lb2+1:N
if strcmp(method2,'joint_sic')
for m=1:M
y_soft_tr(m,k)=filt_seq_tr(m,k);
for m1=1:M
y_soft_tr(m,k)=y_soft_tr(m,k)-symb_hat_soft_tr(m1,k-1:-1:k-Lb2)*B2((m1-1)*(Lb2+1)+2:m1*(Lb2+1),m);
end;
if (m>1)
y_soft_tr(m,k)=y_soft_tr(m,k)-symb_hat_soft_tr(1:m-1,k)'*B2(1:Lb2+1:(m-1)*(Lb2+1),m);
end;
symb_hat_soft_tr(m,k)=tanh(y_soft_tr(m,k)/sqrt(Mse2(m)));
end;
elseif strcmp(method2,'joint_sicrev')
for m=M:-1:1
y_soft_tr(m,k)=filt_seq_tr(m,k);
for m1=1:M
y_soft_tr(m,k)=y_soft_tr(m,k)-symb_hat_soft_tr(m1,k-1:-1:k-Lb2)*B2((m1-1)*(Lb2+1)+2:m1*(Lb2+1),m);
end;
if (m<M)
y_soft_tr(m,k)=y_soft_tr(m,k)-symb_hat_soft_tr(m+1:M,k)'*B2(m*(Lb2+1)+1:(Lb2+1):end,m);
end;
symb_hat_soft_tr(m,k)=tanh(y_soft_tr(m,k)/sqrt(Mse2(m)));
end;
end;
end;

% Diversity combining for hard decision feedback
err_seq=abs(sign(symb-symb_hat));
err_seq_tr=abs(sign(symb_tr-symb_hat_tr));

% Diversity combining for soft decision feedback (BPSK)
err_seq_soft=abs(sign(symb-sign(symb_hat_soft)));
err_seq_soft_tr=abs(sign(symb_tr-sign(symb_hat_soft_tr)));

% Linear combining for hard DFE(Add & detect)
symb_hat_bidfe(:,1:Lb1)=symb(:,1:Lb1);
symb_hat_bidfe(:,N-Lb2+1:N)=symb(:,N-Lb2+1:N);
for m=1:M
symb_hat_bidfe(m,Lb1+1:N-Lb2)=pamdet(Alpha(m)*y(m,Lb1+1:N-Lb2)+(1-Alpha(m))*fliplr(y_tr(m,Lb1+1:N-Lb2)),2,1);
end

% Linear combining for soft DFE(Add & detect)
symb_hat_soft_bidfe(:,1:Lb1)=symb(:,1:Lb1);
symb_hat_soft_bidfe(:,N-Lb2+1:N)=symb(:,N-Lb2+1:N);
for m=1:M
symb_hat_soft_bidfe(m,Lb1+1:N-Lb2)=pamdet(Alpha(m)*y_soft(m,Lb1+1:N-Lb2)+(1-Alpha(m))*fliplr(y_soft_tr(m,Lb1+1:N-Lb2)),2,1);
end;

% Performance measures

ser_ideal(:,iter)=sum(ceil((sign(abs(symb-symb_hat_ideal)-1e-8)+1)/2),2)/size(symb,2);
ser(:,iter)=sum(ceil((sign(abs(symb-symb_hat)-1e-8)+1)/2),2)/size(symb,2);
ser_soft(:,iter)=sum(ceil((sign(abs(symb-sign(symb_hat_soft))-1e-8)+1)/2),2)/size(symb,2);
ser_ideal_tr(:,iter)=sum(ceil((sign(abs(symb_tr-symb_hat_ideal_tr)-1e-8)+1)/2),2)/size(symb_tr,2);
ser_tr(:,iter)=sum(ceil((sign(abs(symb_tr-symb_hat_tr)-1e-8)+1)/2),2)/size(symb_tr,2);
ser_soft_tr(:,iter)=sum(ceil((sign(abs(symb_tr-sign(symb_hat_soft_tr))-1e-8)+1)/2),2)/size(symb_tr,2);
ser_bidfe(:,iter)=sum(ceil((sign(abs(symb-symb_hat_bidfe)-1e-8)+1)/2),2)/size(symb,2);
ser_soft_bidfe(:,iter)=sum(ceil((sign(abs(symb-symb_hat_soft_bidfe)-1e-8)+1)/2),2)/size(symb,2);

end;

toc








