clear;
K=8;
P=179;
L=3;
M=16;
J=11;
U=10^5;

sigma=0.0446    ;

%求出R矩陣書上的11頁(單位圓上取點),通     
R=zeros(M,J);%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

for u=1:M;
    for i=1:J;
        R(u,i)=exp(j*2*pi*((u-1)+(i-1)*M)/(M*J));%*****************************************R為M*J矩陣
    end
end
clear i,j;
R2=R(2,:);                 %                 the second user   1*J 
%//////////////////////////////////////////////////// /////////////

%/////////////////////////////////////////////////////////////
%求出Thita矩陣,書上的12頁,VANDERMONDE MATRIX,通
Thita1=zeros(J,K);

 for i=1:J;
     for j=1:K; 
     Thita1(i,j)=R2(i).^(-j+1);%***************************************************Thita1 為J*K矩陣
     end 
  end
clear i,j;
%////////////////////////////////////////////////////////////////
 

S=randn(K,U);
for i=1:K
    for j=1:U;
        if S(i,j)<0
            S(i,j)=-1;
        else 
            S(i,j)=1;
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%SS1=reshape(S,1,K*U);%????????????????????????????????????????????????????????????????????????????????????

S1=Thita1*S;

%////////////////////////////////////////////////夾在//之間的程序是通的
%求 信道矩陣]
clear i,j;
h=randn(L+1,M)./sqrt(32)+i*randn(L+1,M)./sqrt(32);%********************************一次隨機產生后,應該是一個定值,與信道有關 (L+1)*M,

%///////////////////////////////////////////

%///////////////////////////////////////////////////////////////  only the second user 
Hm=zeros(M,J);%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
for u=1:M;
 for j=1:J;
      %for m=1:M;----------------------------------------------------------原來的錯誤使Hm矩陣出現很多個0
        for l=1:L+1;%????????????????????????????????????????????????????????????//書上是0到L
        Hm1=h(l,u)*exp(R(u,j).^(-l+1));
        Hm(u,j)=Hm(u,j)+Hm1;%*************************************************************************Hm為M*J矩陣
        end
        %end
 end 
end
%////////////////////////////////////////////////////////

%DuHru=zeros(J,J,M);%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

%///////////////////////
%for i=1:M;
%for j=1:J;
%DuHru(j,j,i)=Hm(i,j);%***********************************************************************DuHru為J*J*M矩陣
%end
%end
%DuHru2=DuHru(:,:,2);
%////////////////////////////

for j=1:J; 
DuHru2(j,j)=Hm(2,j);%***********************************************************************DuHru為J*J*M矩陣
 end


S2=DuHru2*S1;

%  NOISE
Gu=zeros(J,P);
 for i=1:J;
     for j=1:P; 
     Gu(i,j)=R(i).^(-j+1);%***************************************************Thita1 為J*K矩陣
     end 
  end
  clear i,j;
 ELTA=sigma.*randn(P,U); 
 %EITA=sigma.*randn(size(Gu));
 noise=Gu*ELTA; 
 S3=S2+noise;
%////////////////////////////////////////////////////////////////


%%%%%%%%%EQUALIZATION]
Geiu1=DuHru2*Thita1;

 Geiu=pinv(Geiu1);
 S4=Geiu*S3;
S5=real(S4);%???????????????????????????????????????????????????????????????????????????????? right?

for i=1:K
    for j=1:U;
        if S5(i,j)<0
            S5(i,j)=-1;
        else 
            S5(i,j)=1;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%COMPARE
biterror=0;
for i=1:K
    for j=1:U;
        if S5(i,j)~=S(i,j)
          biterror=biterror+1;
            
        end
    end
end
biterror
ratio=biterror./(K*U)