% ummse_fse_dfe   Designs Unbiased MMSE T/2-spaced feedforward equalizer (FSE) 
%                 and T-spaced decision feedback equalizer  (DFE) 
%                 for a T/2-spaced complex FIR channel at a certain 
%                 SNR (dB) of AWGN and a unit power iid source given 
%                 a cursor delay and an optional decision feedback delay. 
% 
%                 [f,d,ummse] = ummse_fse_dfe(h,snrdb,Nf,Nd,delta1,delta2)
%   
%                 returns:
%                            f      :   FSE          
%                            d      :   DFE          
%                            ummse  :   Unbiased Minimum Mean Square Error
%                                       (i.e. design is scaled so that
%                                        cursor is set to unity)
%                 given:    
%                            h      :   channel
%                            snrdb  :   SNR in dB
%                            Nf     :   FSE length
%                            Nd     :   DFE length
%                            delta1 :   cursor delay 
%                            delta2 :   decision feedback delay
%                                       (Optional, with default 
%                                        delta2=delta1+1. Need
%                                        delta2>delta1) 


function [f,d,ummse] = ummse_fse_dfe(h,snrdb,Nf,Nd,delta1,delta2)


% Check if delta2 is being used and
% if not, set delta2=delta1+1
if (nargin==5)
   delta2=delta1+1;
elseif (nargin<5)
   error('Missing arguments');
elseif (nargin>6)
   error('Too many arguments');
end;

% Check if delta2>delta1
if (delta2<delta1)
   error('Need delta2>delta1');
end;

% Get SNR from SNR in dB
% and NSR=lambda
snr=10^(snrdb/10);
lambda=1/snr;

% Set up convolution matrix
if (size(h,1)<size(h,2))
   h=h';
end;
HH=convmtx(h,Nf);
H=HH(1:2:size(HH,1),:);

% Channel-FSE combination t=Hf and
% single spike response td
Nt=size(H,1);
td=zeros(Nt,1);
td(delta1)=1;

% Check if delta1 or delta2 are
% two large
if (delta1>Nt) 
   error('delta1 is too large');
elseif (delta2>Nt) 
   error('delta2 is too large');
end;


% Find MMSE-FSE+DFE
I=eye(size(H,2));
Ndmax=min(Nt-delta2+1,Nd);
W=diag([ones(delta2-1,1); zeros(Ndmax,1); ones(Nt-delta2-Ndmax+1,1)]);
f=inv(H'*W*H+lambda*I)*(H'*td);

% Unbias MMSE-FSE+DFE to get UMMSE-FSE+DFE
t=H*f;
f=f/max(t);
t=H*f;
d=[t(delta2:delta2+Ndmax-1); zeros(Nd-Ndmax,1)];
dd=[zeros(delta2-1,1); t(delta2:delta2+Ndmax-1); zeros(Nt-delta2+1-Ndmax,1)];
ummse=norm(t-dd-td)^2+lambda*norm(f)^2;
