adaptive mmse space-time multi-user detection .m

adaptive mmse space-time multiuser detection algorithm.m是自适应MMSE空时多用户检测算法 dec multiuser.m是空时解相关多用户检测

clear;
  usernumber=15; % 小区里的用户数
  Lc1=31; % 用于对符号扩频的长度
  ini(1,:)=[0 0 0 0 1];
  ini(2,:)=[0 0 0 1 0];
  ini(3,:)=[0 0 0 1 1];
  ini(4,:)=[0 0 1 0 0];
  ini(5,:)=[0 0 1 0 1];
  ini(6,:)=[0 0 1 1 0];
  ini(7,:)=[0 0 1 1 1];
  ini(8,:)=[0 1 0 0 0];
  ini(9,:)=[0 1 0 0 1];
  ini(10,:)=[0 1 0 1 0];
  ini(11,:)=[0 1 0 1 1];
  ini(12,:)=[0 1 1 0 0];
  ini(13,:)=[0 1 1 0 1];
  ini(14,:)=[0 1 1 1 0];
  ini(15,:)=[0 1 1 1 1];
  ini(16,:)=[1 0 0 0 0];
  ini(17,:)=[1 0 0 0 1];
  ini(18,:)=[1 0 0 1 0];
  ini(19,:)=[1 0 0 1 1];
  ini(20,:)=[1 0 1 0 0];
  ini(21,:)=[1 0 1 0 1];
  ini(22,:)=[1 0 1 1 0];
  ini(23,:)=[1 0 1 1 1];
  ini(24,:)=[1 1 0 0 0];
  ini(25,:)=[1 1 0 0 1];
  ini(26,:)=[1 1 0 1 0];
  ini(27,:)=[1 1 0 1 1];
  ini(28,:)=[1 1 1 0 0];
  ini(29,:)=[1 1 1 0 1];
  ini(30,:)=[1 1 1 1 0];
  poly1=[1 0 0 1 0 1];
  n=5;
  for k=1:1:usernumber;
      mseq=msequence(ini(k,:),n,poly1);
      PN(k,:)=conver(mseq,n);
  end     % 产生用户PN码
  
   %方向向量的产生
 for i=1:1:usernumber;
     aa1(i)=4*3.14*(rand-0.5);%每个用户第1个路径的波达方向
     aa2(i)=4*3.14*(rand-0.5);%每个用户第2个路径的波达方向
     aa3(i)=4*3.14*(rand-0.5);%每个用户第3个路径的波达方向
     aa4(i)=4*3.14*(rand-0.5);%每个用户第4个路径的波达方向
 end
for i=1:1:usernumber;
    a1(:,i)=[1 exp(j*3.14*sin(aa1(i))) exp(j*2*3.14*sin(aa1(i))) exp(j*3*3.14*sin(aa1(i)))]';%每个用户在天线1上的方向向量
    a2(:,i)=[1 exp(j*3.14*sin(aa2(i))) exp(j*2*3.14*sin(aa2(i))) exp(j*3*3.14*sin(aa2(i)))]';%每个用户在天线2上的方向向量
    a3(:,i)=[1 exp(j*3.14*sin(aa3(i))) exp(j*2*3.14*sin(aa3(i))) exp(j*3*3.14*sin(aa3(i)))]';%每个用户在天线3上的方向向量
    a4(:,i)=[1 exp(j*3.14*sin(aa4(i))) exp(j*2*3.14*sin(aa4(i))) exp(j*3*3.14*sin(aa4(i)))]';%每个用户在天线4上的方向向量
end
  
  
  snrindb=5; %信噪比（分贝） 
  snr=10^(snrindb/10);
%噪声方差假设为1，信号的幅度也设为1
  sgma=1;
%比特能量
  Eb=2*sgma^2*snr;
 %设有15个用户,基站有4个天线，每个用户有4个路径,最大时延为5个码片,同时考虑多址干扰和码片间干扰
  zuidashiyanmapianshu=5;%最大时延码片数
  kuozhanmachang=31+zuidashiyanmapianshu;%扩展码长
 
 %幅度产生
 A(1,1)=0.7*sqrt(Eb);
 A(1,2)=0.2*sqrt(Eb);
 A(1,3)=0.5*sqrt(Eb);
 A(1,4)=0.3*sqrt(Eb);
 for i=2:1:usernumber;
     for j=1:1:4;
         A(i,j)=(2-rand)*sqrt(Eb);;
     end
 end
  
 %各个用户各个路径的时延码片个数
  shiyan1=[3 2 4 5 0 1 5 3 5 2 1 3 1 4 5 5 3 0 0 3 2 1 0 4 3 2 1 4 3 0];%各个用户第1条路径的时延码片数
  shiyan2=[0 3 2 1 5 4 0 2 4 1 2 1 0 5 3 3 5 3 4 0 3 2 4 5 0 1 5 3 5 2];%各个用户第2条路径的时延码片数
  shiyan3=[2 1 0 4 3 2 1 4 3 0 5 4 5 1 0 0 4 4 2 1 4 4 5 5 1 5 4 1 1 5];%各个用户第3条路径的时延码片数
  shiyan4=[4 4 5 5 1 5 4 1 1 5 4 2 2 3 4 1 2 2 3 2 0 3 2 1 5 4 0 2 4 1];%各个用户第4条路径的时延码片数
 
  %扩展码的产生
  for i=1:1:usernumber;
     C1(i,:)=[zeros(1,shiyan1(i)) PN(i,:) zeros(1,zuidashiyanmapianshu-shiyan1(i))];
     if shiyan1(i)==0;
        C1f(i,:) =zeros(1,kuozhanmachang);
    else
       for j=1:1:kuozhanmachang;
           if j<=shiyan1(i);
             C1f(i,j)=PN(i,31-shiyan1(i)+j);
          else C1f(i,j)=0;
          end
       end
     end
     if shiyan1(i)==zuidashiyanmapianshu;
        C1b(i,:) =zeros(1,kuozhanmachang);
     else
       for j=1:1:kuozhanmachang;
           if j<=shiyan1(i)+31;
             C1b(i,j)=0;
          else C1b(i,j)=PN(i,j-31-shiyan1(i));
          end
       end
     end
    
     
     C2(i,:)=[zeros(1,shiyan2(i)) PN(i,:) zeros(1,zuidashiyanmapianshu-shiyan2(i))];
     if shiyan2(i)==0;
        C2f(i,:) =zeros(1,kuozhanmachang);
    else
       for j=1:1:kuozhanmachang;
           if j<=shiyan2(i);
             C2f(i,j)=PN(i,31-shiyan2(i)+j);
          else C2f(i,j)=0;
          end
       end
     end
     if shiyan2(i)==zuidashiyanmapianshu;
        C2b(i,:) =zeros(1,kuozhanmachang);
     else
       for j=1:1:kuozhanmachang;
           if j<=shiyan2(i)+31;
             C2b(i,j)=0;
          else C2b(i,j)=PN(i,j-31-shiyan2(i));
          end
       end
     end
     
     
     C3(i,:)=[zeros(1,shiyan3(i)) PN(i,:) zeros(1,zuidashiyanmapianshu-shiyan3(i))];
     if shiyan3(i)==0;
        C3f(i,:) =zeros(1,kuozhanmachang);
    else
       for j=1:1:kuozhanmachang;
           if j<=shiyan3(i);
             C3f(i,j)=PN(i,31-shiyan3(i)+j);
          else C3f(i,j)=0;
          end
       end
     end
     if shiyan3(i)==zuidashiyanmapianshu;
        C3b(i,:) =zeros(1,kuozhanmachang);
     else
       for j=1:1:kuozhanmachang;
           if j<=shiyan3(i)+31;
             C3b(i,j)=0;
          else C3b(i,j)=PN(i,j-31-shiyan3(i));
          end
       end
     end
     
     
     C4(i,:)=[zeros(1,shiyan4(i)) PN(i,:) zeros(1,zuidashiyanmapianshu-shiyan4(i))];
     if shiyan4(i)==0;
        C4f(i,:) =zeros(1,kuozhanmachang);
    else
       for j=1:1:kuozhanmachang;
           if j<=shiyan4(i);
             C4f(i,j)=PN(i,31-shiyan4(i)+j);
          else C4f(i,j)=0;
          end
       end
     end
     if shiyan4(i)==zuidashiyanmapianshu;
        C4b(i,:) =zeros(1,kuozhanmachang);
     else
       for j=1:1:kuozhanmachang;
           if j<=shiyan4(i)+31;
             C4b(i,j)=0;
          else C4b(i,j)=PN(i,j-31-shiyan4(i));
          end
       end
     end
     
 end

 

 b=zeros(usernumber,10000);   %定义所有用户发送的符号为一个矩阵
 %用户数据的产生 
 for j=1:1:2000+2;
     b(1,j)=1;%用户1的训练序列，假设都为1
 end
 for i=2:1:usernumber;
     for j=1:1:2000+2;
        b(i,j)=sign(0.5-rand);
     end
 end
 
 
%RLS算法产生时间域的权值向量c1、c2、c3、c4,其中c1、c2、c3、c4分别是每个天线在时间域上对应的权值
%初始化
XX=zeros(4,kuozhanmachang,2000);%定义一个矩阵XX
r=0.99;%遗忘因子
P1=0.005*ones(kuozhanmachang,kuozhanmachang)-0.01*rand(kuozhanmachang,kuozhanmachang);%初始化P1
c1=0.005*ones(kuozhanmachang,1)-0.01*rand(kuozhanmachang,1);%初始化权值c1

P2=0.005*ones(kuozhanmachang,kuozhanmachang)-0.01*rand(kuozhanmachang,kuozhanmachang);%初始化P2
c2=0.005*ones(kuozhanmachang,1)-0.01*rand(kuozhanmachang,1);%初始化权值c2

P3=0.005*ones(kuozhanmachang,kuozhanmachang)-0.01*rand(kuozhanmachang,kuozhanmachang);%初始化P3
c3=0.005*ones(kuozhanmachang,1)-0.01*rand(kuozhanmachang,1);%初始化权值c3

P4=0.005*ones(kuozhanmachang,kuozhanmachang)-0.01*rand(kuozhanmachang,kuozhanmachang);%初始化P4
c4=0.005*ones(kuozhanmachang,1)-0.01*rand(kuozhanmachang,1);%初始化权值c4
for l=2:1:2000;
    % 产生每个天线的时间域接收信号
    X=zeros(4,kuozhanmachang);
    for i=1:1:usernumber;
        X=X+a1(:,i)*A(i,1)*C1(i,:)*b(i,l)+a1(:,i)*A(i,1)*C1b(i,:)*b(i,l+1)+a1(:,i)*A(i,1)*C1f(i,:)*b(i,l-1)+a2(:,i)*A(i,2)*C2(i,:)*b(i,l)+a2(:,i)*A(i,2)*C2b(i,:)*b(i,l+1)+a2(:,i)*A(i,2)*C2f(i,:)*b(i,l-1)+a3(:,i)*A(i,3)*C3(i,:)*b(i,l)+a3(:,i)*A(i,3)*C3b(i,:)*b(i,l+1)+a3(:,i)*A(i,3)*C3f(i,:)*b(i,l-1)+a4(:,i)*A(i,4)*C4(i,:)*b(i,l)+a4(:,i)*A(i,4)*C4b(i,:)*b(i,l+1)+a4(:,i)*A(i,4)*C4f(i,:)*b(i,l-1);   
    end
    X=X+sgma*randn(4,kuozhanmachang);
    XX(:,:,l)=X;%将接收数据存储到XX中
    %天线1时间域上的RLS算法
    K1=(1/(r+X(1,:)*P1*X(1,:)'))*P1*X(1,:)';
    e1=1-X(1,:)*c1;
    c1=c1+K1*e1;
    P1=(1/r)*(P1-K1*X(1,:)*P1);
    %天线2时间域上的RLS算法
    K2=(1/(r+X(2,:)*P2*X(2,:)'))*P2*X(2,:)';
    e2=1-X(2,:)*c2;
    c2=c2+K2*e2;
    P2=(1/r)*(P2-K2*X(2,:)*P2);
    %天线3时间域上的RLS算法
    K3=(1/(r+X(3,:)*P3*X(3,:)'))*P3*X(3,:)';
    e3=1-X(3,:)*c3;
    c3=c3+K3*e3;
    P3=(1/r)*(P3-K3*X(3,:)*P3);
    %天线4时间域上的RLS算法
    K4=(1/(r+X(4,:)*P4*X(4,:)'))*P4*X(4,:)';
    e4=1-X(4,:)*c4;
    c4=c4+K4*e4;
    P4=(1/r)*(P4-K4*X(4,:)*P4);
end

%RLS算法产生空间域的权值向量w1
%初始化
r=0.99;%遗忘因子
P=0.005*ones(4,4)-0.01*rand(4,4);%初始化P
w1=0.005*ones(4,1)-0.01*rand(4,1);%初始化w1
for l=2:1:2000;
    % 产生每个天线的接收信号
    X=XX(:,:,l);%接收数据
    Y1=X(1,:)*c1;%天线1时域滤波输出信号
    Y2=X(2,:)*c2;%天线2时域滤波输出信号
    Y3=X(3,:)*c3;%天线3时域滤波输出信号
    Y4=X(4,:)*c4;%天线4时域滤波输出信号
    Y=[Y1 Y2 Y3 Y4]';%将Y1、Y2、Y3、Y4组成一个向量
    %空域中的RLS算法
    K=(1/(r+Y'*P*Y))*P*Y;
    e=1-Y'*w1;
    w1=w1+K*e;
    P=(1/r)*(P-K*Y'*P);
end



%重新产生用户数据
for i=1:1:usernumber;
     for j=1:1:10000+2;
        b(i,j)=sign(0.5-rand);
     end
 end

%接收信号并计算误码率
error=0;%误码数置0
for l=2:1:10000;
    % 产生每个天线的时间域接收信号
    X=zeros(4,kuozhanmachang);
    for i=1:1:usernumber;
        X=X+a1(:,i)*A(i,1)*C1(i,:)*b(i,l)+a1(:,i)*A(i,1)*C1b(i,:)*b(i,l+1)+a1(:,i)*A(i,1)*C1f(i,:)*b(i,l-1)+a2(:,i)*A(i,2)*C2(i,:)*b(i,l)+a2(:,i)*A(i,2)*C2b(i,:)*b(i,l+1)+a2(:,i)*A(i,2)*C2f(i,:)*b(i,l-1)+a3(:,i)*A(i,3)*C3(i,:)*b(i,l)+a3(:,i)*A(i,3)*C3b(i,:)*b(i,l+1)+a3(:,i)*A(i,3)*C3f(i,:)*b(i,l-1)+a4(:,i)*A(i,4)*C4(i,:)*b(i,l)+a4(:,i)*A(i,4)*C4b(i,:)*b(i,l+1)+a4(:,i)*A(i,4)*C4f(i,:)*b(i,l-1);   
    end
    X=X+sgma*randn(4,kuozhanmachang);
    %时域滤波
    Y1=X(1,:)*c1;
    Y2=X(2,:)*c2;
    Y3=X(3,:)*c3;
    Y4=X(4,:)*c4;
    %空域滤波
    Y=[Y1 Y2 Y3 Y4]';
    YY=w1'*Y;
    b1=sign(real(YY));%取YY的实部，并得到对用户1数据的估计
    if b1~=b(1,l);
        error=error+1;%当出现误码时，error加1
    end
end
wumalv=error/10000;%计算误码率