% DFEperformance.m 
% Shengli Zhou, 11/20/2002

f1=[0.04, 0.05, 0.07, 0.21, 0.5, 0.72, 0.36, 0.21, 0.03, 0.06]';
f2=[0.407, 0.815, 0.407]';
f3=[0.227, 0.460, 0.688, 0.460, 0.227]';

EbN0dB=[0:2:35];
Eb=1;
N0s=Eb./(10.^(EbN0dB/10));
K1=15;	 
% K2 is assumed to be equal to L

for chan=1:3
    if chan ==1
       f=f1;
    elseif chan ==2
      f=f2;
    elseif chan ==3
      f=f3;
    end
    for isnr=1:length(EbN0dB);
   	N0=N0s(isnr);
	L=length(f)-1;
	Psi=zeros(K1+1,K1+1);
	for ix=1:K1+1
	    for iy=1:K1+1
	    l=ix-K1-1;
	    j=iy-K1-1;
	    tp=0;
	    for m=0:(-1*l)
                if m <= L & m+l-j >=0 & m+l-j <=L 
                   tp=tp+conj(f(m+1))*f(m+l-j+1);
		end
	    end
	    Psi(ix,iy)=tp;
	    end
	end
        % BPSK constellation, we restrict the noise to be real
	Psi=Psi+N0/2*eye(K1+1);
        xi=[zeros(K1-L,1); flipud(f)]; 
        rho2=xi'*(Psi\xi);
	SNR=rho2/(1-rho2);
	BER(chan,isnr)=Qfunc(sqrt(SNR));
	BER_awgn(isnr)=Qfunc(sqrt(2*1/N0));
    end
end
        
figure(1)
semilogy(EbN0dB,BER(3,:),'k-o'); hold on
semilogy(EbN0dB,BER(2,:),'k-s'); hold on
semilogy(EbN0dB,BER(1,:),'k--'); hold on
semilogy(EbN0dB,BER_awgn,'k'); hold on
hold off
xlabel('SNR');
ylabel('Probability of Error');
axis([0 35 1e-4 4e-1]);
legend('Channel C','Channel B','Channel A','AWGN');