function [wzfe_ff, wzfe_fb,sigma2_zfedfe, SINR_zfedfe, wmmse_ff,wmmse_fb, sigma2_mmsedfe, SINR_mmsedfe]=dfe_equalizer(p,L,Lb,delta, sigma2_noise)
% [wzfe_ff, wzfe_fb,sigma2_zfedfe, SINR_zfedfe, wmmse_ff,wmmse_fb,sigma2_mmsedfe, SINR_mmsedfe]=dfe_equalizer(p,L,Lb,delta, sigma2_noise)
% ZFE and MMSE decision feedback equalizer function
% Inputs
%  p - sampled channel pulse response vector
%  L - feedforward linear equalizer length
%  Lb - feedback linear equalizer length
%  delta - delay of the feedforward equalizer
%  sigma2_noise - noise variance
%
% Outputs
%  wzfe_ff - feedforward zero-forcing filter coefficients
%  wzfe_fb - feedback zero-forcing filter coefficients
%  sigma2_zfedfe - mean square error ZFE-DFE
%  SINR_zfedfe - Signal-to-Interference and Noise Ratio for ZFE-DFE
%  wmmse_ff - feedforward mmse filter coefficients
%  wmmse_fb - feedback mmse filter coefficients
%  sigma2_mmsedfe - mean square error MMSE-DFE
%  SINR_mmsedfe - Signal-to-Interference and Noise Ratio for MMSE-DFE

% Make channel Toeplitz matrix
P=convmtx(p',L)';

% Make a delay target vector
delta_vec=[zeros(1,delta) 1 zeros(1,size(P,2)-delta-1)];

Pdfe=[P(:,1:delta+1) P(:,delta+1+Lb+1:end)];
delta_vec_ff=[delta_vec(1:delta+1) delta_vec(delta+1+Lb+1:end)];   
wzfe_ff_check=delta_vec_ff*Pdfe'*(Pdfe*Pdfe'+1e-12*eye(size(Pdfe,1)))^-1;
wmmse_ff_check=delta_vec_ff*Pdfe'*(Pdfe*Pdfe'+sigma2_noise*eye(size(Pdfe,1)))^-1;
% alternatively you can define 
Jd=[zeros(1,delta+1) ones(1,Lb) zeros(1,size(P,2)-1-delta-Lb)]';

wzfe_ff=delta_vec*P'*(P*P'-P*Jd*Jd'*P')^-1;
wmmse_ff=delta_vec*P'*(P*P'-P*Jd*Jd'*P'+sigma2_noise*eye(size(P,1)))^-1;
wzfe_fb=wzfe_ff*P*Jd;
wmmse_fb=wmmse_ff*P*Jd;

sigma2_zfedfe=1-wzfe_ff*P*delta_vec'+sigma2_noise*sum(wzfe_ff.^2);
sigma2_mmsedfe=1-wmmse_ff*P*delta_vec';
SINR_zfedfe=(wzfe_ff*P*delta_vec')^2/sigma2_zfedfe;
SINR_mmsedfe=1/sigma2_mmsedfe-1;