subspacemud.m

此程序为DS-CDMA链路中加入基于子空间的PASTd多用户检测方法。系统中采用m序列阔频

function [ich10,qch10,ERRORo,ERRORs]=Subspacemud(ich8,qch8,MUDmode,code)
switch MUDmode
    case 1
        ich10=ich8;
        qch10=qch8;
    case 2
        ich10=ich8;
        qch10=qch8;
    case 3
        ich10=ich8;
        qch10=qch8;
    case 4
 %********************************PASTd*****************************       
        [l,m]=size(code);
        [n,v]=size(ich8);
        ndata= v/m;
        X=zeros(m,n+1);   %initialazing parameters 
        y=zeros(1,n);
        ich10=zeros(1,ndata);
        namuda=zeros(1,n);
        namudapre=rand(1,n);
        Us=zeros(m,n);
%         Uspretemp=ones(1,m)*5;
%         Uspre=diag(Uspretemp);
        Uspre=rand(m,n);
        beta = 1;
        ich8temp=ich8';
        qch8temp=qch8';
        cv=ones(1,m);
        C=diag(cv);
        ERRORo=zeros(1,ndata);
        ERRORs=zeros(1,ndata);
        for i = 1:ndata
            X(:,1)=ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1);
            for k = 1:n             %begin to estimation
%             fid=fopen('PASTdU.mat','r');
%             [Upre,count]=fread(fid,[v,n],int8);
%             fclose(fid);
%             fid=fopen('PASTdna.mat','r')
%             [namudapre,count]=fread(fid,n,int8);
%             fclose(fid);
            y(k)=Uspre(:,k)'*X(:,k);
            namuda(k)=beta*namudapre(k)+abs(y(k))^2;
            Us(:,k)=Uspre(:,k)+(X(:,k)-Uspre(:,k)*y(k))*conj(y(k))/namuda(k);
            X(:,k+1)=X(:,k)-Us(:,k)*y(k);
            end
            yest=Us.'*(ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1));
            codetemp=code';
            p1est=Us.'*codetemp(:,1);
            c=zeros(1,n);
            for ii=1:n
            c(ii)=1/namuda(1,ii);
            end
            c1=diag(c)*p1est;
            ich10(i)=real(c1'*yest);
            qch10(i)=imag(c1'*yest);
            Uspre=Us;
            namudapre=namuda;
            
%%%%%%%%%%%%%%%%%测试正交错误
            vtesto=ones(1,n);   
            Etesto=diag(vtesto);
            Errtesto=Etesto-Us'*Us;
            ErrFtesto=0;
            for iii=1:n
                for jjj=1:n
                   ErrFtesto=ErrFtesto+abs(Errtesto(iii,jjj)^2);
               end
            end
            ERRORo(i)=20*log10(sqrt(ErrFtesto)/sqrt(n));
            
%%%%%%%%%%%%%%%%%测试子空间错误            
            C=C+(ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1))*(ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1))';
            [Px,D]=eigs(C,n);
            vtests=ones(1,m);
            Etests=diag(vtests);
            Errtests=(Etests-Us*inv(Us'*Us)*Us')*(Px*Px');
            ErrFtests=0;
            for iii=1:m
                for jjj=1:m
                   ErrFtests=ErrFtests+abs(Errtests(iii,jjj)^2);
               end
            end
            ERRORs(i)=20*log10(sqrt(ErrFtests)/sqrt(n));
        end
 
  %*********************************************PAST*************************      
    case 5 
        [l,m]=size(code);
        [n,v]=size(ich8);
        ndata= v/m;
        W=zeros(m,n);
        Wpretemp=ones(1,m);
        Wpretemp1=diag(Wpretemp);
        Wpre=Wpretemp1(:,1:n);
%         Wpre=rand(m,n);
        Zpretemp=ones(1,n)
        Zpre=diag(Zpretemp);
        ich8temp=ich8';
        qch8temp=qch8';
        alpha=1;
        beta=1;
        namudaPre=zeros(1,n);
        cv=ones(1,7);
        C=diag(cv);
        ERRORo=zeros(1,ndata);
        ERRORs=zeros(1,ndata);
        for i =1:ndata
           r=ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1);
           y=Wpre'*r;
           q=Zpre*y;                       %h
%            gama=1/(1+y'*q);               %%%%%%%%%%%%%%%%
%            p=gama*(r-Wpre*y);
%            W=Wpre+p*q';                   %OPAST文档中的方法
%            Z=Zpre-gama*q*q'               %%%%%%%%%%%%%%%%
           g=q/(1+y'*q);                    %g%%%%%%%%%%%%%%%%% 
           Z=Zpre-g*q';
           p=r-Wpre*y;                      %e PAST的方法 
           W=Wpre+p*g';                     %%%%%%%%%%%%%%%%%%
           ytemp=y.';
           for ii=1:n
              namuda(ii)=beta*namudaPre(ii)+abs(ytemp(ii))^2;
           end
           Wpre=W;
           namudaPre=namuda;
           yest=W.'*(ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1));
           codetemp=code';
           p1est=W.'*codetemp(:,1);
           c=zeros(1,n);
           for ii=1:n
           c(ii)=1/namuda(1,ii);
           end
           c1=diag(c)*p1est;
           ich10(i)=real(c1'*yest);
           qch10(i)=imag(c1'*yest);
           
%%%%%%%%%%%%%%%%%测试正交错误
            vtesto=ones(1,n);   
            Etesto=diag(vtesto);
            Errtesto=Etesto-W'*W;
            ErrFtesto=0;
            for iii=1:n
                for jjj=1:n
                   ErrFtesto=ErrFtesto+abs(Errtesto(iii,jjj)^2);
               end
            end
            ERRORo(i)=20*log10(sqrt(ErrFtesto)/sqrt(n));
            
            
            %%%%%%%%%%%%%%%%%测试子空间错误            
            C=C+(ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1))*(ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1))';
            [Px,D]=eigs(C,n);
            vtests=ones(1,m);
            Etests=diag(vtests);
            Errtests=(Etests-W*inv(W'*W)*W')*(Px*Px');
            ErrFtests=0;
            for iii=1:m
                for jjj=1:m
                   ErrFtests=ErrFtests+abs(Errtests(iii,jjj)^2);
               end
            end
            ERRORs(i)=20*log10(sqrt(ErrFtests)/sqrt(n));
            
            
       end


%*********************************************OPAST*************************      
    case 6 
        [l,m]=size(code);
        [n,v]=size(ich8);
        ndata= v/m;
        W=zeros(m,n);
        Wpretemp=ones(1,m)*0.1;
        Wpretemp1=diag(Wpretemp);
        Wpre=Wpretemp1(:,1:n);
%         Wpre=rand(m,n);
        Zpretemp=ones(1,n)*3;
        Zpre=diag(Zpretemp);
        ich8temp=ich8';
        qch8temp=qch8';
        alpha=1;
        beta=1;
        namudaPre=zeros(1,n)
        cv=ones(1,7);
        C=diag(cv);
        ERRORo=zeros(1,ndata);
        ERRORs=zeros(1,ndata);
        for i =1:ndata
           r=ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1);
           y=Wpre'*r;
           q=(1/alpha)*Zpre*y;
           gama=1/(1+y'*q);
           p=gama*(r-Wpre*y);
           sita=(1/(q'*q))*(1/(sqrt(1+(p'*p)*(q'*q)))-1);
           pO=sita*Wpre*q+(1+sita*(q'*q))*p;
           W=Wpre+pO*q';
           Z=(1/alpha)*Zpre-gama*q*q';
           ytemp=y.';
           for ii=1:n
              namuda(ii)=beta*namudaPre(ii)+abs(ytemp(ii))^2;
           end
           Wpre=W;
           Zpre=Z;
           namudaPre=namuda;
           yest=W.'*(ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1));
           codetemp=code';
           p1est=W.'*codetemp(:,1);
           c=zeros(1,n);
           for ii=1:n
           c(ii)=1/namuda(1,ii);
           end
           c1=diag(c)*p1est;
           ich10(i)=real(c1'*yest);
           qch10(i)=imag(c1'*yest);
           
%%%%%%%%%%%%%%%%%测试正交错误
            vtesto=ones(1,n);   
            Etesto=diag(vtesto);
            Errtesto=Etesto-W'*W;
            ErrFtesto=0;
            for iii=1:n
                for jjj=1:n
                   ErrFtesto=ErrFtesto+abs(Errtesto(iii,jjj)^2);
               end
            end
            ERRORo(i)=20*log10(sqrt(ErrFtesto)/sqrt(n));
            
   %%%%%%%%%%%%%%%%%测试子空间错误            
            C=C+(ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1))*(ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1))';
            [Px,D]=eigs(C,n);
            vtests=ones(1,m);
            Etests=diag(vtests);
            Errtests=(Etests-W*inv(W'*W)*W')*(Px*Px');
            ErrFtests=0;
            for iii=1:m
                for jjj=1:m
                   ErrFtests=ErrFtests+abs(Errtests(iii,jjj)^2);
               end
            end
            ERRORs(i)=20*log10(sqrt(ErrFtests)/sqrt(n));
       end
       %********************************OPASTd*****************************       
   case 7 
        [l,m]=size(code);
        [n,v]=size(ich8);
        ndata= v/m;
        X=zeros(m,n+1);   %initialazing parameters 
        y=zeros(1,n);
        ich10=zeros(1,ndata);
        namuda=zeros(1,n);
        namudapre=rand(1,n);
        Us=zeros(m,n);
%         Uspretemp=ones(1,m)*5;
%         Uspre=diag(Uspretemp);
        Uspre=rand(m,n);
        beta = 1;
        ich8temp=ich8';
        qch8temp=qch8';
        Wcov=zeros(1,ndata);
        btemp=zeros(m,1);
        b=zeros(m,n);
        cv=ones(1,7);
        C=diag(cv);
        ERRORo=zeros(1,ndata);
        ERRORs=zeros(1,ndata);
        for i = 1:ndata
            X(:,1)=ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1);
            for k = 1:n             %begin to estimation
%             fid=fopen('PASTdU.mat','r');
%             [Upre,count]=fread(fid,[v,n],int8);
%             fclose(fid);
%             fid=fopen('PASTdna.mat','r')
%             [namudapre,count]=fread(fid,n,int8);
%             fclose(fid);
            y(k)=Uspre(:,k)'*X(:,k);
            namuda(k)=beta*namudapre(k)+abs(y(k))^2;
            Us(:,k)=Uspre(:,k)+(X(:,k)-Uspre(:,k)*y(k))*conj(y(k))/namuda(k);
%             if i<200              %%%%%%%%%%%%次数规定
               if k==1              %%%%%%%%%%%%%%%正交化%%%%%%%%%%%
                  b(:,k)=Us(:,k);
               else
                 for kk=1:k-1
                    btemp=btemp+((b(:,kk)'*Us(:,k))/(b(:,kk)'*b(:,kk)))*b(:,kk);
                 end
                 b(:,k)=Us(:,k)-btemp;
                 btemp=zeros(m,1);
               end                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
               X(:,k+1)=X(:,k)-b(:,k)*y(k);
%              else                          %%%%%%%%%次数规定
%                X(:,k+1)=X(:,k)-Us(:,k)*y(k);
%              end                          %%%%%%%%%次数规定
            end                       
%             for k=1:n          %%%%%%%%%%%%%%%%%%%%单位化%%%%%%%%%%%%%%%%
%                 b(:,k)=b(:,k)/sqrt(b(:,k)'*b(:,k));
%             end                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            Us=b;
            yest=Us.'*(ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1));
            codetemp=code';
            p1est=Us.'*codetemp(:,1);
            c=zeros(1,n);
            for ii=1:n
            c(ii)=1/namuda(1,ii);
            end
            c1=diag(c)*p1est;
            ich10(i)=real(c1'*yest);
            qch10(i)=imag(c1'*yest);
            Uspre=Us;
            namudapre=namuda;
            
            %%%%%%%%%%%%%%%%%测试正交错误
            vtesto=ones(1,n);   
            Etesto=diag(vtesto);
            Errtesto=Etesto-Us'*Us;
            ErrFtesto=0;
            for iii=1:n
                for jjj=1:n
                   ErrFtesto=ErrFtesto+abs(Errtesto(iii,jjj)^2);
               end
            end
            ERRORo(i)=20*log10(sqrt(ErrFtesto)/sqrt(n));
            
%%%%%%%%%%%%%%%%%测试子空间错误            
            C=C+(ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1))*(ich8temp((i-1)*7+1:i*7,1)+j*qch8temp((i-1)*7+1:i*7,1))';
            [Px,D]=eigs(C,n);
            vtests=ones(1,m);
            Etests=diag(vtests);
            Errtests=(Etests-Us*inv(Us'*Us)*Us')*(Px*Px');
            ErrFtests=0;
            for iii=1:m
                for jjj=1:m
                   ErrFtests=ErrFtests+abs(Errtests(iii,jjj)^2);
               end
            end
            ERRORs(i)=20*log10(sqrt(ErrFtests)/sqrt(n));
        
            
     
        end
 end