clear
fix(clock)
dim=10;
  N=29; %粒子數(shù)
iter=30; %迭代次數(shù)
R3=(1/31)*[31 -1 -9 7 -1 7 7 -9 -1 -1
    -1 31 7 -9 -1 7 -9 7 -1 -1
    -9 7 31 -1 7 -1 -1 -1 7 -9
    7 -9 -1 31 7 -1 -1 -1 7 7
    -1 -1 7 7 31 7 -9 -9 -1 -1
    7 7 -1 -1 7 31 -1 -1 -9 7
    7 -9 -1 -1 -9 -1 31 -1 7 7
    -9 7 -1 -1 -9 -1 -1 31 7 7
    -1 -1 7 7 -1 -9 7 7 31 -1
    -1 -1 -9 7 -1 7 7 7 -1 31]; %10個用戶

Numpso=zeros(1,7);Numomd=zeros(1,7);Numcd=zeros(1,7);Numde=zeros(1,7);


De=1:1024; %for de:=1 to 16 產(chǎn)生 此段程序是用來產(chǎn)生所有可能的輸入來進行最優(yōu)檢測
inpu=zeros(10,1024);	%產(chǎn)生10×1024的零矩陣
for d1=1:1024		
   DEm=De(1,d1);		%按d1順序取de中的值
   for d2=1:10			
     inpu(d2,d1)=rem(DEm,2);	%求余mod(dem,2)
     DEm=fix(DEm/2);  % 返回整數(shù)部分
   end
end

for d3=1:1024
   for d4=1:10
      Trans=inpu(d4,d3);
      if Trans==0
         inpu(d4,d3)=-1;
      end
   end
end


Rate=[1 1.6 2.5 4 6.3 10 15.849];%0,2,4,6,8,10,12dB


for    t=1:7
    
    W=sqrt(eye(dim)*Rate(t));  
    H=zeros(dim,dim);
    H=W'*R3*W;       
    input=zeros(dim,N);
    for i=1:N%第一個是傳統(tǒng)檢測器輸出值    
        for j=1:dim
            input(j,i)=rand(1)-0.5;   
        end
    end
    
    for i=1:30   %1000
        b0=sign(rand(1,dim)-0.5);
        if b0==0
            b0=-1
        end
        b0=b0';
        %產(chǎn)生噪聲n（10，1）方差為1
        n=randn(1,dim);
        n=n';
        %產(chǎn)生輸出Y
         Y=H*b0+n;          %%Y=R*W*b0+n;

        Bcd=sign(Y) ;%此段程序是進行CD檢測即傳統(tǒng)檢測器檢測
        Numcd(t)=Numcd(t)+sum(ne(Bcd,b0));

        R4=inv(R3);  %此段程序是進行去相關檢測
        Bde=sign(R4*Y);%(10,1)
        Numde(t)=Numde(t)+sum(ne(Bde,b0));

           z=[];
          for i=1:2.^10;
         z=[z Y'*inpu(:,i)-0.5*inpu(:,i)'*H*inpu(:,i)];
        end;
        [vm I]=max(z);
        Bomd=inpu(:,I);% 最佳多用戶檢測結果
        Numomd(t)=Numomd(t)+sum(ne(Bomd,b0));



 %此段程序是關于基本粒子群優(yōu)化算法
        E1=zeros(1,N);
        v1=zeros(dim,N);
        c1=zeros(dim,N);
        c1=input;
        v1=2*0.1*rand(dim,N)-0.1;  %%v1=rand(10,N)-0.5;
        for l=1:iter %終止條件迭代次數(shù)30    
            for i1=1:N
                b1=c1(:,i1);         
                E1(i1)=0.5*b1'*H*b1-b1'*Y; % 評價各個粒子的函數(shù)值 E1(i1)=0.5*b1'*H*b1-Y'*b1;
            end
            %初始局部最小函數(shù)值即單個粒子所達到的最小函數(shù)值p及其對應位置bp
            if l==1    
                p1=zeros(1,N);
                bp1=zeros(dim,N);
                p1=E1;
                bp1=c1;     
                [g1 i2]=min(E1);%初始全局最小函數(shù)值g及其所對應的位置bg
                bg1=zeros(dim,1);             
                bg1=c1(:,i2); 
            end
            for i3=1:N
                if E1(i3)<p1(i3)
                    p1(i3)=E1(i3);
                    bp1(:,i3)=c1(:,i3);
                end
            end
            for i4=1:N
                if p1(i4)<g1
                    g1=p1(i4);
                    bg1=bp1(:,i4);               %m=i4;
                end
            end
            %根據(jù)下式改變每個粒子的速度與位置
            for i5=1:N
                
                b1=c1(:,i5);
                v1(:,i5)=v1(:,i5)+2*rand(dim,1).*(bp1(:,i5)-b1)+ 2*rand(dim,1).*(bg1-b1);
                for i6=1:dim
                    for j2=1:N
                        if v1(i6,j2)<-0.1  %-0.5
                            v1(i6,j2)=-0.1; %-0.5
                        end
                        if v1(i6,j2)>0.1  %0.5
                            v1(i6,j2)=0.1;  %0.5
                        end
                    end
                end
                b1=b1+v1(:,i5);
                if b1<-0.5
                   b1=-0.5;
                end
                if b1>0.5
                   b1=0.5;
                end
                c1(:,i5)=b1;
            end       
        end
        pso=sign(bg1);
        Numpso(t)=Numpso(t)+sum(ne(pso,b0));
                

                 end
            end
            Numcd=Numcd/300, 
            Numde=Numde/300,
            Numpso=Numpso/300,
            Numomd= Numomd/300,
     x=[0 2 4 6 8 10 12];
        plot(x,Numcd,'b-*',x,Numde,'r-o',x,Numpso,'m-x',x,Numomd,'k:'),
         
    
       


            fix(clock)