% Bidirectional DFE
% With MMSE-BiDFE optimized DFE taps

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


% Generate Channel & equalizer

C1=[0.3162 1 0.3162];
C=C1;
C=C./sqrt(C*C');
C_est=C; %Estimate of channel to be used for resconstruction

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

[f_opt,b_opt,delay,mmse]=trdfe_opt_taps(C,Lf,Lb,10^(-snr(snr_index)/10));
% use last input as 0 for unconstrained tap optimization

Lc=length(C);
Lbad=1;%choice of delay for BAD comparison, a fn of largest tap of C

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);

% 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

% Feedforward filtering
rec_seq_tr=fliplr(rec_seq);
rec_seq1_tr=conv(rec_seq_tr,transpose(f_opt));
r_tr=rec_seq1_tr(delay: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. 

%symb_hat_up=zeros(1,P*N-P+1);
%symb_hat_up(1:P:P*N-P+1)=symb_hat;

r_hat_ideal=conv(symb_hat_ideal,C_est);
r_hat_ideal_tr=conv(fliplr(symb_hat_ideal_tr),C_est);
r_hat=conv(symb_hat,C_est);
r_hat_tr=conv(fliplr(symb_hat_tr),C_est);
%r_hat_soft=conv(symb_hat_soft,C_est);
%r_hat_soft_tr=conv(fliplr(symb_hat_soft_tr),C_est);

r_diff_ideal=(r_hat_ideal-rec_seq).^2;
r_diff_ideal_tr=(r_hat_ideal_tr-rec_seq).^2;
r_diff=(r_hat-rec_seq).^2;
r_diff_tr=(r_hat_tr-rec_seq).^2;
%r_diff_soft=(r_hat_soft-rec_seq).^2;
%r_diff_soft_tr=(r_hat_soft_tr-rec_seq).^2;

symb_bad_ideal(1:N_avg)=symb_hat_ideal(1:N_avg);
symb_bad_ideal(N:-1:N_avg+1)=symb_hat_ideal_tr(1:N-N_avg);
symb_bad(1:N_avg)=symb_hat(1:N_avg);
symb_bad(N:-1:N_avg+1)=symb_hat_tr(1:N-N_avg);
%symb_bad_soft(1:N_avg)=symb_hat_soft(1:N_avg);
%symb_bad_soft(N:-1:N_avg+1)=symb_hat_soft_tr(1:N-N_avg);

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

if (sum(r_diff_ideal(k-N_avg:1:k+N_avg)) < sum(r_diff_ideal_tr(k-N_avg:1:k+N_avg))) 
symb_bad_ideal(k-Lbad)=symb_hat_ideal(k-Lbad);
end;

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;

%if (sum(r_diff_soft(k-N_avg:1:k+N_avg)) < sum(r_diff_soft_tr(k-N_avg:1:k+N_avg))) 
%symb_bad_soft(k-Lbad)=symb_hat_soft(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((y(Lb+1:N-Lb)+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((y_soft(Lb+1:N-Lb)+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_ideal(iter)=symberror(symb,symb_bad_ideal);
ser_bad(iter)=symberror(symb,symb_bad);
%ser_bad_soft(iter)=symberror(symb,sign(symb_bad_soft));

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_ideal(snr_index)=mean(ser_bad_ideal)
SER_bad(snr_index)=mean(ser_bad)
%SER_bad_soft(snr_index)=mean(ser_bad_soft)

end;

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

figure(1);
pl1=semilogy(snr,SER_ideal,'r');
hold on;
pl2=semilogy(snr,SER,'b-.');
pl3=semilogy(snr,SER_bidfe,'k-^');
pl4=semilogy(snr,SER_bidfe_soft,'g-*');
xlabel('SNR  in  dB','FontSize',18);
ylabel('Symbol  Error  Rate','FontSize',18);
legend([pl1,pl2,pl3,pl4],'Ideal feedback','Normal Mode DFE','LC-BiDFE: hard','LC-BiDFE: soft');
grid;

toc

save results_c1_opt 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 SER_bad_ideal; 





