% Bidirectional Arbitrated DFE 
% Implements the BAD algorithm proposed by McGahey et al.

% Initialization
rand('seed',0);
block_num=1000*ones(1,20);   
% Number of blocks
N=1000;           % block length
Lf=8;            % Number of feedfroward taps
Lb=4;             % Number of feedback taps
P=1;              % Fractional sampling factor
snr=[4:1:16];            % in dB
N_avg=5;        % moving average window length on each side for BAD

% Generate Channel & equalizer

C1=[0.3162 1 0.3162];
C2=[0.2 1 0.3162 -0.2 0.1 -0.05 0.02];
C3=[-0.3920  0.0152  0.8320  0.0152  -0.3920]

C=C2;
C=C./sqrt(C*C');
Lc=length(C);

Lbad=1;%choice of delay for BAD comparison, a fn of largest tap of C
C_est=C;

tic

for snr_index=1:length(snr)
snr(snr_index)
sig=10^(-snr(snr_index)/20);

% Channel estimation and DFE tap estimation
[f,b,delay,mmse]=dfe_tap_estimator(C,Lf,Lb,10^(-snr(snr_index)/10),P,0);
[f_tr,b_tr,delay_tr,mmse_tr]=dfe_tap_estimator(fliplr(C),Lf,Lb,10^(-snr(snr_index)/10),P,0);
[mmse_div,rho,alpha]=trdfe_mse(C_est,f,b,delay,f_tr,b_tr,delay_tr,10^(-snr(snr_index)/10));

for iter=1:block_num(snr_index)

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

% Generate Source Sequence
symb=pam(N,2,1); 
symb_up=zeros(1,P*N-P+1);
symb_up(1:P:P*N-P+1)=symb;

% Transmit through channel
rec_seq=conv(symb_up,C);
rec_seq=rec_seq+noise(length(rec_seq),snr(snr_index),1);

f_opt=f;b_opt=b;

% Feedforward filtering
rec_seq1=conv(rec_seq,transpose(f_opt));
r=rec_seq1(delay:P:end); % downsampled signal

% Feedback filtering

% Ideal feedback
symb_hat_ideal(1:Lb)=symb(1:Lb);
for k=Lb+1:N
y(k)=r(k)-symb(k-1:-1:k-Lb)*b_opt;
symb_hat_ideal(k)=pamdet(y(k),2,1);
end;

% Hard Decision feedback
symb_hat(1:Lb)=symb(1:Lb);
for k=Lb+1:N
y(k)=r(k)-symb_hat(k-1:-1:k-Lb)*b_opt;
symb_hat(k)=pamdet(y(k),2,1);
end;

% Soft Decision feedback
symb_hat_soft(1:Lb)=symb(1:Lb);
for k=Lb+1:N
y_soft(k)=r(k)-symb_hat_soft(k-1:-1:k-Lb)*b_opt;
symb_hat_soft(k)=tanh(y_soft(k)/sig);
end;

% Time reversed decision feedback
f_opt=f_tr;b_opt=b_tr;

% Feedforward filtering
rec_seq_tr=fliplr(rec_seq);
rec_seq1_tr=conv(rec_seq_tr,transpose(f_opt));
r_tr=rec_seq1_tr(delay_tr:P:end); % downsampled signal

% Feedback filtering 

% Ideal feedback
symb_tr=fliplr(symb);
symb_hat_ideal_tr(1:Lb)=symb_tr(1:Lb);
for k=Lb+1:N
y_tr(k)=r_tr(k)-symb_tr(k-1:-1:k-Lb)*b_opt;
symb_hat_ideal_tr(k)=pamdet(y_tr(k),2,1);
end;

% Hard Decision feedback
symb_hat_tr(1:Lb)=symb_tr(1:Lb);
for k=Lb+1:N
y_tr(k)=r_tr(k)-symb_hat_tr(k-1:-1:k-Lb)*b_opt;
symb_hat_tr(k)=pamdet(y_tr(k),2,1);
end;

% Soft Decision feedback
symb_hat_soft_tr(1:Lb)=symb_tr(1:Lb);
for k=Lb+1:N
y_soft_tr(k)=r_tr(k)-symb_hat_soft_tr(k-1:-1:k-Lb)*b_opt;
symb_hat_soft_tr(k)=tanh(y_soft_tr(k)/sig);
end;

%%%%% Reconstruction of received Sequence
% Note: r_hat should be compared with rec_seq and not r.
% Currently, will work only for baud-spaced. 

r_hat=conv(symb_hat,C_est);
r_hat_tr=conv(fliplr(symb_hat_tr),C_est);
r_diff=(r_hat-rec_seq).^2;
r_diff_tr=(r_hat_tr-rec_seq).^2;

symb_bad(1:N_avg)=symb_hat(1:N_avg);
symb_bad(N:-1:N_avg+1)=symb_hat_tr(1:N-N_avg);

for k=N_avg+1+Lbad:N-N_avg+Lbad,

if (sum(r_diff(k-N_avg:1:k+N_avg)) < sum(r_diff_tr(k-N_avg:1:k+N_avg))) 
symb_bad(k-Lbad)=symb_hat(k-Lbad);
end;

end;

%%%%%%%%End Reconstruction

%%%%%%%%%%Begin BiDFE

% Suboptimal combining for hard DFE(Add & detect)
symb_bidfe(1:Lb)=symb(1:Lb);
symb_bidfe(N-Lb+1:N)=symb(N-Lb+1:N);
symb_bidfe(Lb+1:N-Lb)=pamdet((alpha*y(Lb+1:N-Lb)+(1-alpha)*fliplr(y_tr(Lb+1:N-Lb)))/2,2,1);

% Suboptimal combining for soft DFE(Add & detect)
symb_bidfe_soft(1:Lb)=symb(1:Lb);
symb_bidfe_soft(N-Lb+1:N)=symb(N-Lb+1:N);
symb_bidfe_soft(Lb+1:N-Lb)=pamdet((alpha*y_soft(Lb+1:N-Lb)+(1-alpha)*fliplr(y_soft_tr(Lb+1:N-Lb)))/2,2,1);


% Plots & Performance measures

ser_ideal(iter)=symberror(symb,symb_hat_ideal);
ser(iter)=symberror(symb,symb_hat);
ser_soft(iter)=symberror(symb,sign(symb_hat_soft));
ser_ideal_tr(iter)=symberror(symb_tr,symb_hat_ideal_tr);
ser_tr(iter)=symberror(symb_tr,symb_hat_tr);
ser_soft_tr(iter)=symberror(symb_tr,sign(symb_hat_soft_tr));
ser_bidfe(iter)=symberror(symb,symb_bidfe);
ser_bidfe_soft(iter)=symberror(symb,symb_bidfe_soft);
ser_bad(iter)=symberror(symb,symb_bad);

end;

SER_ideal(snr_index)=mean(ser_ideal)
SER(snr_index)=mean(ser)
SER_soft(snr_index)=mean(ser_soft)
SER_ideal_tr(snr_index)=mean(ser_ideal_tr)
SER_tr(snr_index)=mean(ser_tr)
SER_soft_tr(snr_index)=mean(ser_soft_tr)
SER_bidfe(snr_index)=mean(ser_bidfe)
SER_bidfe_soft(snr_index)=mean(ser_bidfe_soft)
SER_bad(snr_index)=mean(ser_bad)

end;

set(0,'DefaultAxesFontSize',18);
set(0,'DefaultTextFontSize',18);
set(0,'DefaultTextFontName','times');
set(0,'DefaultAxesFontName','times');

figure(1);
pl1=semilogy(snr,SER,'b-.','LineWidth',2);
hold on;
pl2=semilogy(snr,SER_tr,'c-v','LineWidth',2);
pl3=semilogy(snr,SER_bidfe,'k-o','LineWidth',2);
pl4=semilogy(snr,SER_bidfe_soft,'g-*','LineWidth',2);
pl5=semilogy(snr,SER_bad,'m-^','LineWidth',2);
xlabel('SNR  in  dB','FontSize',18);
ylabel('Symbol  Error  Rate','FontSize',18);
legend([pl1,pl2,pl3,pl4,pl5],'Normal mode DFE','Time-Reversal Mode DFE','LC-BiDFE (Hard)','LC-BiDFE (Soft)','BAD');
grid;

toc

save results_ser_bad_c5 block_num snr N_avg N C Lf Lb P Lbad SER_ideal SER SER_soft SER_ideal_tr SER_tr SER_soft_tr SER_bidfe SER_bidfe_soft SER_bad; 





