gmc.m

简单实用的编码GMC的曲线

%主程序
clear;
 clf;
clc;
L=3;
M=16;
J=10;
K=8;
P=143;

%////////////////////////////////////////////////夹在//之间的程序是通的
%求 信道矩阵
h=ones(L+1,M)/0.60;%**********************************************************************一次随机产生后,应该是一个定值,与信道有关,
%///////////////////////////////////////////




%求出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
%/////////////////////////////////////////////////////////////////


%/////////////////////////////////////////////////////////////
%求出Thita矩阵,书上的12页,VANDERMONDE MATRIX,通
Thita=zeros(J,K,M);

for k=1:M;
 for i=1:J;
     for j=1:K; 
     Thita(i,j,k)=R(k,i).^(j-1);%***************************************************Thita 为J*K*M矩阵
     end 
  end
end
%////////////////////////////////////////////////////////////////



%///////////////////////////////////////////////////////////////
Hm=zeros(M,J);%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
for u=1:M;  
    for i=1:J;
     Hm(u,i)=0;
     end
end

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));
        Hm(u,j)=Hm(u,j)+Hm1;%*************************************************************************Hm为M*J矩阵
        end
        %end
 end 
end
%////////////////////////////////////////////////////////

 
%////////////////////////////////////////////////////////////
%第u个用户的第k个码的能量,通
for u=1:M;
    for k=1:K;
        E(u,k)=M*J;
    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
%////////////////////////////

 
%////////////////////////////////////////////////////////////////////// 
Geiu1=zeros(J,K,M);

for i=1:M;
geiu=DuHru(:,:,i)*Thita(:,:,i);
 Geiu1(:,:,i)=geiu;%*************************************************************************全部问题的关键
end
%/////////////////////////////////////////////////////////////// 
    
 
%/////////////////////////////////
Geiu=zeros(K,J,M);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

for i=1:M;
Geiu(:,:,i)=pinv(Geiu1(:,:,i));%伪逆    ************************************************************为K*J*M矩阵
end
%//////////////////////////////////这个循环可插入上一个循环中

   
%Geiu已知 ///////////////////////////////////////////////////////

Gu=zeros(J,P,M);
for i=1:M;
    for k=1:K;
        for j=1:J;
Gu(j,k,i)=Thita(j,k,i);%**************************************************************************Gu为J*P*M矩阵
end
end
end
%
%



GGU=zeros(K,P,M);%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

for i=1:M;
GGU(:,:,i)=Geiu(:,:,i)*Gu(:,:,i);% ********************************************************GGU为K*P*M矩阵
end
%////////////////////////////////////////////////////////



%///////////////////////////////////////////////////////////
%以下为画图测试,这里用nn取代Eb/No,通
nn=10.001:0.5:14.001;
Eb=1;
No=1;
Pe2=0;
for u=1:M;%
    for k=1:K;
      Pe1= Qfunct( real(sqrt(1/(GGU(k,:,u)*(GGU(k,:,u)').*E(u,k)))*sqrt(2*10.^(nn/10))));%real后面的其实是实数,但是是x+0i的形式,所以取实数根本没变化值
      Pe2=Pe2+Pe1;
    end
end
Pe=Pe2/(M*J);

    Pe4=0;
 for u=1:M;%
    for k=1:K;
      Pe3= (2^11-1)*Qfunct( 11/15*3*real(sqrt(1/(GGU(k,:,u)*(GGU(k,:,u)').*E(u,k)))*sqrt(2*10.^(nn/10))));%real后面的其实是实数,但是是x+0i的形式,所以取实数根本没变化值
      Pe4=Pe3+Pe4;
    end
end
Pe5=Pe4/(M*J);




semilogy(nn,Pe,'ko:');
hold on
semilogy(nn,Pe5,'k*:');
grid on;
legend('M=16,K=8,L=3','coded GMC-CDMA',3);
title(['GMC-CDMA' ]);
xlabel('Eb/No[dB]');
ylabel('average BER');
hold off;
%/////////////////////////////////////////////////////////////