function [wzfe,sigma2_zfe, sigma2_zfe_res_isi, SINRzfe, wmmse, sigma2_mmse, sigma2_mmse_res_isi, SINRmmse]=linear_equalizer(p,L,delta, sigma2_noise)
% [wzfe,sigma2_zfe, sigma2_zfe_res_isi, SINRzfe, wmmse, sigma2_mmse, sigma2_mmse_res_isi, SINRmmse]=linear_equalizer(p,L,delta, sigma2_noise)
% ZFE and MMSE linear equalizer function
% Inputs
%  p - sampled channel pulse response vector
%  L - feedforward linear equalizer length
%  delta - delay of the feedforward equalizer
%  sigma2_noise - noise variance
%
% Outputs
%  wzfe - zero-forcing filter coefficients
%  sigma2_zfe - mean square error ZFE
%  sigma2_zfe_res_isi - residual ISI ZFE error
%  SINRzfe - Signal-to-Interference and Noise Ratio for ZFE
%  wmmse -  mmse filter coefficients
%  sigma2_mmse - mean square error MMSE
%  SINRmmse - Signal-to-Interference and Noise Ratio for MMSE

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

% Zero forcing solution
wzfe=delta_vec*P'*(P*P')^-1;

% calcualte equalized pulse response
peqzfe=wzfe*P;

% calculate the mean-square error
sigma2_zfe_res_isi = 1-wzfe*P*delta_vec';
sigma2_zfe_noise_enh = sigma2_noise*sum(wzfe.^2);
sigma2_zfe = sigma2_zfe_res_isi + sigma2_zfe_noise_enh;
SINRzfe=max(peqzfe)^2/sigma2_zfe; % in books it says Ex/sigma2_zfe, but our ZFE is finite length, hence cannot perfectly equalize the signal

% figure(1)
% hold off;
% subplot(211), stem(wzfe);
% subplot(212), stem(peqzfe);
% disp(wzfe);
% disp(peqzfe);

%%% Do the MMSE linear equalizer
wmmse=delta_vec*P'*(P*P'+sigma2_noise*eye(size(P,1)))^-1;

% calcualte equalized pulse response
peqmmse=wmmse*P;

% calculate the mean-square error
sigma2_mmse= 1-wmmse*P*delta_vec';
sigma2_mmse_res_isi=sum(peqmmse.^2)-max(peqmmse)^2;
SINRmmse_check=1/sigma2_mmse-1; % unbiased SINRmmse
SINRmmse=max(peqmmse)^2/(sigma2_mmse_res_isi+sigma2_noise*sum(wmmse.^2)); %%% this definition is more intuitive since we did not specifically cover biased vs. unbiased SINR
%%% note that the SINR expressions above are equivalent



% figure(2)
% hold off;
% subplot(211), stem(wmmse);
% subplot(212), stem(peqmmse);
% disp(wmmse);
% disp(peqmmse);