% MMSEperformance.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));
K=15;
Length=2*K+1;

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;
        x=xcorr(f); x=x(:)';
        Row=[x(L+1:end), zeros(1,Length-L-1)];
        % BPSK constellation, we restrict the noise to be real
        Gam=toeplitz(Row)+N0/2*eye(Length);
        xi=[zeros(K-L,1); flipud(f); zeros(K,1)];
        rho2=xi'*(Gam\xi);
	SNR=rho2/(1-rho2);
	BER(chan,isnr)=Q(sqrt(SNR));
	BER_awgn(isnr)=Q(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');