N=64;K=32;z=[6 7 8 9 10 11 12 13 14 15 16 17];snr=power(10,z/10);R=[1 2 3 4 5 6 7 8];L=power(2,R);delta=[1.5956,0.9957,.5860,0.3352,0.1881,0.1041,0.0569,0.0308];for i=1:N    for j=1:K        r=randint;        s(i,j)=(2*r-1)/sqrt(N);    endendfor i=1:8    x(1)=-inf;    x(L(i)+1)=inf;            for k=2:L(i)        x(k)=(k-(L(i)+2)/2)*delta(i);    end    for k=1:L(i)         y(k)=(k-(L(i)+1)/2)*delta(i);    end    g=0;    for k=1:L(i)         g=g+y(k)*(exp(-x(k)^2/2)-exp(-x(k+1)^2/2));    end    alpha_g(i)=g/(sqrt(2*pi));    e=0;    for k=1:L(i)        e=e+(y(k)^2)*(qfunc(x(k))-qfunc(x(k+1)));    end    var_eta(i)=e-alpha_g(i)^2;end        a=s(1:N,5);    for h=1:8    for j=1:12     c=1+var_eta(h)*(1+K*snr(j)/N)/(alpha_g(h)^2);                b=inv(snr(j)*s*eye(K,K)*s'+eye(N,N)*c);        U(h,j)=snr(j)*a'*b*a;        gamma(h,j)=U(h,j)/(1-U(h,j));    ber(h,j)=qfunc(sqrt(gamma(h,j)));    end    end%     hold on    for h=1:5      ber1(h,:)= log10(ber(h,1:12));%        plot(z,ber1)    endplot(z,ber1(1,:),'--r');grid on;hold on;plot(z,ber1(2,:),'--b');grid on;hold on;plot(z,ber1(3,:),'--m');grid on;hold on;plot(z,ber1(4,:),'--k');grid on;hold on;plot(z,ber1(5,:),'--g');grid on;hold on;