% function []= eye_eq(h,Ex,noise_var,eq_type);
%-----------------------------------------------------------------------------
%  INPUTS
% h         = pulse response, example  [a*D^-1 + 1] ==> [a 1]
% Ex        = average energy of signals, Ex_bar
% noise_var = noise variance (sigma squared)
% eq_type   = Z => ZERO FORCING (may have numerical issues, so use M with zero noise instead
%             M => MMSE
%             D => MMSE-DFE
% OUTPUTS, some are graphs
% outputs: pe_(zfe/mmle_dfe/mmse_le)= probability of error with Equalizer from N input 2PAM[+1/-1] data sequences  
% outputs: pe_no_eq  = probability of error without Equalizer
% outputs: dfseSNR   = receiver(equalizer)SNR, unbiased in dB
% outputs: pe_SNR = error probability estimation from receiver(equalizer)SNR, Pe = Q function of sqrt(dfseSNR)
% this function shows Frequency response of the channel p and equalizer filter, eye diagram, receiver SNR, and probability of error.
% N = 250; % # of input 2 PAM[+1/-1] data, you can increase the N of input data sequences in order to calculate the Pe accurately
% created 1/06 by Jungsub Byun and M. Malkin  EE379A
% -------------------------------------------------------------------------
function []= eye_eq(h,Ex,noise_var,eq_type);
close all; 
warning off all
fltrLength = 40;  %default feedforward equalizer length 40
size = length(h);
delay = -1;  %search for optimum delay
nff = max(fltrLength,3*length(h));
eq_type = lower(eq_type);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N =  250; % # of input 2 PAM[+1/-1] data, you can increase the N of input data sequences in order to calculate the Pe accurately 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nbb = 0;

if eq_type == 'z'  %zero-forcing
    for i=1:size,
        p_z(i)=h(size-i+1); % (D = Z^-1), 1+a*D^-1=1+a*Z ==>[1 a] 
    end
    if abs(roots(p_z)) > 1
        error('channel has zeros outside unit circle and can not be inverted with ZF Equalizer');
    else
        W_0 = impz(1,p_z); % zf equalizer filter from Z transform,
        nff = length(W_0);
        delay=0;
        dfseSNR = 10*log10(Ex/(noise_var*sum(abs(W_0).^2)));
    end
elseif eq_type == 'm' %mmse-le dfseSNR, W_0
    [dfseSNR, W_0,delay] = dfecolor_m(1,h,nff,nbb,delay,Ex,noise_var*[1 zeros(1,nff-1)]); % mmse-le equalizer filter and SNR
elseif eq_type == 'd' %mmse-dfe dfseSNR, W_0
    nbb = max(2*length(h),floor(fltrLength/2));
    [dfseSNR, W_0,delay] = dfecolor_m(1,h,nff,nbb,delay,Ex,noise_var*[1 zeros(1,nff-1)]); % mmse-sfe equalizer filter and SNR
else
    error('Equalizer type must Z(ZF), M(MMSE), or D(MMSE-DFE)');
end

sigma=sqrt(noise_var);	% standard deviation of the noise
for i=1:N,				% generate input data sequence
    if (rand>0.5),
        I(i)=1;         % Decision is "+1".
    else
        I(i)=-1;        % Decision is "-1".
    end;
end;

%channel p:
A=1;                            % dummy for filter(p,A,I)
rec_sig=filter(h,A,I);    		% the received signal without noise for channel p, A=1, I=generated input data sequence
%% since we want eye diagram scale on [-1,1] divide the noise by
%% sqrt(signal's energy)
noise_gauss = 1/(sqrt(Ex))*sigma*randn(1,N);
y=rec_sig+noise_gauss;	   	    % received signal with noise for channel p

if eq_type == 'z'              %zero-forcing
    W=W_0(1:nff);               % feedforward filter nff taps
    y_zfe=conv(W,y);           % the received signal after ZFE feedforward filter nff taps  
else
    W=W_0(1:nff);            % feedforward filter nff taps
    b=1*W_0(nff+1:nff+nbb);  % feedback filter nbb taps
    
    y_dfe=filter(W,A,y);    % the received signal after feedforward filter nff taps
    
    % detection, and probability of error calculation 
    % after feedforward filter nff taps and feedback filter nbb taps
    numoferr1 = 0;
    r=zeros(1,N);
    decis1=zeros(1,N);
    y_decis_dfe=zeros(1,size+3);
    for i=delay+1:N-1
        y_dfe(delay+1:delay+5)=I(1:5);  % to initialize and compensate for the delay in filtering 
        decis1(delay+1:delay+5)=I(1:5); % to initialize and compensate for the delay in filtering 
        if (y_dfe(i)<0),
            decis1(i)=-1;		      		% Decision is "-1".
        else
            decis1(i)=1;		      		% Decision is "1".
        end;
        decis1_temp=zeros(1,nbb);
        k=i;
        for j=1:nbb
            decis1_temp(j)=decis1(k);
            k=k-1;
        end;
        y_dfe(i+1)=y_dfe(i+1)+ dot(b,decis1_temp); % the operation part of feedback filter nbb taps 
        
        y_decis_dfe(i)=decis1(i);
        if not(y_decis_dfe(i)==I(i-delay))      	% If it is an error, increase the error counter.
            numoferr1=numoferr1+1;                  % accumulating error
        end;
    end;
    
end

%plot frequency response of the channel p
figure;
freqz(h); title('Channel Frequency Response');

%%PLOT EQUALIZATION FILTER RESPONSE & AND EQUALIZED RESPONSE
if eq_type == 'z' % zfe
    figure;freqz(W_0); % ZFE frequency reponse
    title('ZF EQUALIZER FREQUENCY RESPONSE');
    if imag(W_0)==0
        figure;plot(W_0); % ZFE impulse response
        title('ZF EQUALIZER IMPULSE RESPONSE VS # of TAPS');
    end
    figure; freqz(1);
    title(strcat('ZF EQUALIZED CHANNEL ,  SNR = ', num2str(dfseSNR),' dB'));
    numoferr1 = 0;
    r=zeros(1,N);
    decis1=zeros(1,N);
    y_decis_zfe=zeros(1,50);
    for i=delay+1:N-1
        if (y_zfe(i)<0),
            decis1(i)=-1;		      		% Decision is "-1".
        else
            decis1(i)=1;		      		% Decision is "1".
        end;
        y_decis_zfe(i)=decis1(i);
        if not(y_decis_zfe(i)==I(i-delay))      	% If it is an error, increase the error counter.
            numoferr1=numoferr1+1;                  % accumulating error
        end;
    end;
    pe_zfe=numoferr1/(N-delay)  % probability of error estimate after zfe equalizing 
    
elseif eq_type == 'm' %mmse-le
    figure;freqz(W);% mmse-le frequency reponse
    title('MMSE EQUALIZER FREQUENCY RESPONSE');
    if imag(W_0)==0
        figure;plot(W_0);% mmse-le impulse reponse
        title('MMSE EQUALIZER IMPULSE RESPONSE VS # of TAPS');
    end
    figure; freqz(conv(W,h));
    title(strcat('MMSE EQUALIZED CHANNEL , SNR = ', num2str(dfseSNR),' dB'));
    pe_mmse_le=numoferr1/(N-delay) % probability of error estimate after equalizing - feedforward filter mmse-le  
elseif eq_type == 'd' %mmse-dfe
    figure;freqz(W);% mmse-dfe frequency reponse
    title('FEEDFORWARD SECTION MMSE-DFE EQUALIZER FREQUENCY RESPONSE');
    if imag(W_0)==0
        figure;plot(W_0);% mmse-dfe impulse reponse
        title('FEEDFORWARD SECTION MMSE-DFE EQUALIZER IMPULSE RESPONSE VS # of TAPS');
    end
    b_inv = impz(1,b)';
    figure;freqz(conv([zeros(1,delay) b_inv],conv(W_0,h)));
    title(strcat('MMSE-DFE EQUALIZED CHANNEL, SNR = ', num2str(dfseSNR),'dB'));
    pe_mmse_dfe=numoferr1/(N-delay) % probability of error estimate after equalizing - feedforward filter nff taps and feedback filter nbb taps
end
%probability of error estimate without dfe equalizer
numoferr2=0;
decis2=zeros(1,N);
y_decis=zeros(1,N+size+1);
for i=3:N-1,
    if (y(i)<0),
        decis2(i)=-1;		      		% Decision is "-1".
    else
        decis2(i)=1;		      		% Decision is "1".
    end;
    y_decis(i)=decis2(i);
    if (y_decis(i)~=I(i)),    		    % If it is an error, increase the error counter.
        numoferr2=numoferr2+1;          % accumulating error
    end;
end;
pe_no_eq=numoferr2/(N-4)	  	      		% probability of error estimate without dfe equalizer

dfseSNR                                     % outputs: dfseSNR = equalizer SNR, unbiased in dB
if imag(sqrt(dfseSNR))==0
    pe_SNR=q(sqrt(dfseSNR))    % probability of error estimate from dfseSNR, Pe = Q function of sqrt(dfseSNR)
end

% EYE Diagram- before Equalizer (Channel output)
figure;
title('EYE Diagram- before Equalizer (Channel output)')
h=y;
x2=linspace(0,2,2000);
a1=(sinc(x2));
a2=(sinc(x2-1));
a3=(sinc(x2-2));
hold on
for i=delay+1:N-3
    plot(x2, h(i)*a1+h(i+1)*a2+h(i+2)*a3)
    axis([0 2 -1.3 1.3])
end
hold off
% EYE Diagram- After Equalizer
if eq_type == 'z'
    figure;
    title('EYE Diagram - After ZFE')
    h=y_zfe;
elseif eq_type == 'm'
    figure;
    title('EYE Diagram - After MMSE-LE')
    h=y_dfe;
elseif eq_type == 'd'
    figure;
    title('EYE Diagram - After MMSE-DFE ')
    h=y_dfe;
end
x2=linspace(0,2,2000);
a1=(sinc(x2));
a2=(sinc(x2-1));
a3=(sinc(x2-2));
hold on
for i=delay+1:N-3
    plot(x2, h(i)*a1+h(i+1)*a2+h(i+2)*a3)
    axis([0 2 -1.3 1.3])
end
hold off    
warning on all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% DEFINE THIS NEW SUBFUNCTION
%%% ITS EXACTLY LIKE STANDARD EQUALIZE DESIGN PROGRAM
%%% EXCEPTS IT RETURNS OPTIMUM DELAY AS WELL.
%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dfseSNR,w_t,opt_delay]=dfecolor_m(l,p,nff,nbb,delay,Ex,noise);
% ------------------------------------------------------
% l     = oversampling factor
% p     = pulse response, oversampled at l (size)
% nff   = number of feedforward taps
% nbb   = number of feedback taps
% delay = delay of system <= nff+length of p - 2 - nbb
%         if delay = -1, then choose best delay
% Ex    = average energy of signals
% noise = noise autocorrelation vector (size l*nff)
% NOTE: noise is assumed to be stationary
%
% outputs:
% dfseSNR = equalizer SNR, unbiased in dB
% w_t = equalizer coefficients [ w b]
% opt_delay = optimal delay found/delay given
% created 4/96 ; Corrected Dimitrios Bisias, 3/2005; 
% modified 2/06 M. Malkin;
% ------------------------------------------------------

size = length(p);
nu = ceil(size/l)-1;
p = [p zeros(1,(nu+1)*l-size)];

% error check
if nff<=0
    error('number of feedforward taps > 0');
end
if delay > (nff+nu-1-nbb)
    error('delay must be <=  (nff+(length of p)-2-nbb)');
elseif delay < -1 
    error('delay must be >= 0');
elseif delay == -1
    delay = 0:1:nff+nu-1-nbb;
end

%form ptmp = [p_0 p_1 ... p_nu] where p_i=[p(i*l) p(i*l-1)... p((i-1)*l+1)
ptmp(1:l,1) = [p(1); zeros(l-1,1)];
for i=1:nu
    ptmp(1:l,i+1) = conj((p(i*l+1:-1:(i-1)*l+2))');
end

%form matrix P, vector channel matrix
P = zeros(nff*l+nbb,nff+nu);
for i=1:nff,
    P(((i-1)*l+1):(i*l),i:(i+nu)) = ptmp;
end


%precompute Rn matrix - constant for all delays
Rn  = zeros(nff*l+nbb);
Rn(1:nff*l,1:nff*l) = toeplitz(noise);

dfseSNR = -100;
P_init = P;
%loop over all possible delays
for d = delay,
    P = P_init;
    P(nff*l+1:nff*l+nbb,d+2:d+1+nbb) = eye(nbb);
    %P
    temp= zeros(1,nff+nu);
    temp(d+1)=1;
    %construct matrices
    Ry  = Ex*P*P' + Rn;
    Rxy = Ex*temp*P';
    new_w_t = Rxy*inv(Ry);
    sigma_dfse = Ex - new_w_t*Rxy';
    new_dfseSNR = 10*log10(Ex/sigma_dfse - 1);
    %save setting of this delay if best performance so far
    if new_dfseSNR >= dfseSNR
        w_t = new_w_t;
        dfseSNR = new_dfseSNR;
        opt_delay = d;        
    end
end

function [y]=Q(x)
% [y]=Q(x)
%		Q evaluates the Q-function. 
%   		y = 1/sqrt(2*pi) * integral from x to inf of exp(-t^2/2) dt.
%     		y = (1/2) * erfc(x/sqrt(2)).
y=(1/2)*erfc(x/sqrt(2));