clear
echo off

% Adapt a FSE+DFE combination for square M-QAM over 
% T/2-spaced microwave channels with AWGN initialized by  
% Unbiased MMSE-FSE receiver (i.e. MMSE receiver scaled to
% set cursor to unity) 

% Load channel
load /home/combine/raulc/research/sims/comms/channels/chan3.mat
h=C/norm(C);
Nh=length(h);

% Select FSE and DFE lengths
Nf=96;
Nd=32;

% Add nu*T delay for decision feedback
% to channel-FSE combination delay
nu=2;

% Simulation length
L=500000;

% Set M-QAM constellation size
M=64;
a=qam(M,L);
au=zeros(2*L,1);  % Upsample
au(1:2:2*L)=a;

% Define noise & SNR simulation length
snrdb=35;
snr=10^(snrdb/10);
n=sqrt(1/(2*snr))*(randn(2*L,1)+i*randn(2*L,1));

% Transmit
r=filter(h,1,au)+n;

% Initialize FSE to MMSE-FSE and DFE to zero
[f,delta,umse] = ummse_fse(h,snrdb,Nf);
d=zeros(Nd,1);

% Step-size
mu=3e-3;

% Initialize auxiliary variables
ahat=zeros(L,1);
e=zeros(L,1);
z=zeros(L,1);

% Store in memory every K iterations 
K=1000;
fhist=zeros(Nf,L/K);
dhist=zeros(Nd,L/K);


% Adapt
for k=(Nf+Nh):L,

   % FSE and DFE regressors
   R=r(2*k-1:-1:2*k-Nf);
   Ahat=ahat(k-1-nu:-1:k-Nd-nu);
  
   % Find soft and hard estimates 
   z(k)=R.'*f-Ahat.'*d;
   ahat(k)=unqam(z(k),M);

   % Update
   e(k)=z(k)-ahat(k);
   f=f-mu*conj(R)*e(k);
   d=d+mu*conj(Ahat)*e(k);

   % Store parameter trajectories 
   if (rem(k,K)==0)
      fhist(:,k/K)=f;
      dhist(:,k/K)=d;
      disp(k);
   end;

end;

% Find channel-FSE combo
tt=conv(h,f);
t=tt(1:2:length(tt));
Nt=length(t);

% Find effective system response and achieve MSE 
tau=zeros(Nt,1);
tau(delta)=1;
dd=[zeros(delta+nu,1); d; zeros(Nt-Nd-delta-nu,1)];
mse=norm(t-tau-dd)^2+1/snr*norm(f)

% Plot results
% Parameter trajectories
figure(1)
subplot(211);
plot(abs(fhist'));
xlabel('k (x1000)'); ylabel('|f_i(k)|'); 
title('FSE (magnitude) parameter trajectories'); 
subplot(212);
plot(abs(dhist'));
xlabel('k (x1000)'); ylabel('|d_i(k)|'); 
title('DFE (magnitude) parameter trajectories'); 


% Responses
figure(2)
subplot(311);
stem(real(f)); title('(b)');
subplot(312);
stem(real(t)); title('(c)');
subplot(313);
stem(real(dd)); axis([0 200 -0.2 1]);
title('(d)');
orient tall
