
% [w, b, SNRmfb, SNRdmt_teq, SNRdmt, Nstar]=teq(h,L,nu,N,noise,Ex_bar,gap,cb)
%
% programmed by J.M. Cioffi 
% Calculates a time domain equalizer and calculates the SNR for an infinite
% length DMT system to compare
% 
% Inputs
% ------
% h: channel impulse response
% (FIR channels: p= 1+0.5D+2.3D^2->p=[1 0.5 2.3] 
% IIR channels:  p=1/(2.4-0.5D)-> p=[2.4 -0.5])
% L: number of equalizer taps (should be < N/2)
% nu: target channel pulse response length < size(h)-1
% N:  FFT size for eventual computation of SNRdmt_teq
% noise: noise variance per sample
% Ex_bar: Input energy per sample
% cb = 2 for real, 1 for complex
%
% Outputs
% -------
% w: TEQ equalizer
% b: MMSE target channel shaping function (length=nu+1)
% SNRmfb: matched filter bound SNR
% SNRdmt_teq: DMT SNR for this equalized channel for N input (dB scale)
% SNRdmt: DMT SNR for the unequalized first nu+1 samples of the channel
%         Nbig=10000 (dB scale)
% Nopt:  Number of used dimensions
%
% Remarks
% -------
% calls "DMTra.m" and "waterfill_gn" functions
% 
% *****************************************************


function [w, b, SNRmfb, SNRdmt_teq, SNRdmt, Nopt]=teq(h,L,nu,N,noise,Ex_bar,gap,cb)
%   initialization
    norm_h=norm(h);
    [~ , size_h]=size(h);
    if nu > size_h - 2
        'error - nu < size of h - 1'
        return
    elseif L > N/2-1
        'error L > N/2-1'
    else
    C=[h(1) ; zeros(L-1,1)];
    R=[h,zeros(1,L-1)];
    H=toeplitz(C,R);
    btemp=zeros(L,nu+1);
    wtemp=zeros(L,L);
    mmsetemp=zeros(1,L);
    
%   find the MMSE's for all delays
    for delay = 1:L
    ryy=Ex_bar*H*H'+noise*eye(L,L);
    rxy=Ex_bar*[zeros(nu+1,delay) eye(nu+1) zeros(nu+1,L+size_h-2-nu-delay)]*H';
    rle=eye(nu+1)*Ex_bar-rxy*inv(ryy)*rxy';
    [v d]=eig(rle);
    [mmsetemp(delay) , index] = min(diag(d));
    btemp(delay,:)=norm_h*v(:,index)';
    wtemp(delay,:)=btemp(delay,:)*rxy*inv(ryy);
    mmsetemp(delay)=d(index,index)*norm_h^2;
    end
    
%   extract TEQ for this best delay    
    [mmse_opt , index_opt] = min(mmsetemp);
    b=btemp(index_opt,:);
    w=wtemp(index_opt,:);
    eqh=conv(w,h);
    
% compute performance indications
    norm_c=norm(eqh(index_opt+1:index_opt+nu+1));
    error=[eqh(1:index_opt),zeros(1,nu+1),eqh(index_opt+nu+2:end)];
    Re=conv(error,conj(error(end:-1:1)));
    isi_energy=norm(error)^2;
    noise_energy=norm(w)^2*noise;
    Rw=noise*conv(w,conj(w(end:-1:1)));
    distortion=isi_energy*Ex_bar+noise_energy;
    Rb=conv(b,conj(b(end:-1:1)));

    SNRmfb=norm_h^2*Ex_bar/noise;
    SNRmfb=10*log10(SNRmfb);
    ReF=abs(fft(Re,N));
    RwF=abs(fft(Rw,N));
    RbF=abs(fft(Rb,N));
    gn=RbF./(ReF+RwF);
    Nbig=10000;
    fixfac=(Nbig/(Nbig + nu))*((N+nu)/N);
    [bn, en, Nopt] = waterfill_gn(gn,Ex_bar, gap, cb);
    bsumbar=sum(bn)/(N+nu);
    SNRdmt_teq=10*log10(2^(cb*bsumbar)-1);

    %SNRdmt for an infinite length DMT without the equalizer but taking the
    %first nu+1 samples of the channel
    [gn,en_bar,bn_bar,Nstar,b_bar,SNRdmt]=DMTra(h, noise, fixfac*Ex_bar,10000,gap);
 
    end


