Nt=4;%發(fā)射天線數(shù)為4
    Nr=1;%接收天線數(shù)為1
    min_snr=0;
    max_snr=0;%最大的信噪設(shè)定為15
    symbol_per_frame=120;%每一幀的符號數(shù)為120個
    tot_frame_tx=zeros(1,(max_snr-min_snr)/2+1);
    tot_err_symbol =zeros(1,(max_snr-min_snr)/2+1);
    tot_err_bit=zeros(1,(max_snr-min_snr)/2+1);
   for snr=min_snr:2:max_snr
       display(snr);
        for ii=1:10
             tot_frame_tx((snr-min_snr)/2+1) = tot_frame_tx((snr-min_snr)/2+1)+1;
             bit_sequence = rand(1,symbol_per_frame * 2)>0.5;        %生成偽隨機二進制待發(fā)送序列
             mapping_sequence = Mapping(bit_sequence); %二進制信道與mapping信號的轉(zhuǎn)換
            sig =sqrt(1/(10^(snr/10)));  %歸一化的噪聲能量       
            addtional_gaussian_noise =sig * (randn(Nr,Nt,symbol_per_frame+1)/sqrt(2) + j*randn(Nr,Nt,symbol_per_frame+1)/sqrt(2)); %生成加性高斯噪聲
           H=RayleighCH(Nr,Nt); %生成具有瑞利分布特性的信道衰落系數(shù)
            X1=eye(4);
            R1=H*X1;%+addtional_gaussian_noise(:,:,symbol_per_frame+1);      %參考信號  
            for k=1:4:symbol_per_frame-3
               X2=[mapping_sequence(k) -conj(mapping_sequence(k+2)) -conj(mapping_sequence(k+3)) mapping_sequence(k+1); 
                   mapping_sequence(k+2) conj(mapping_sequence(k)) -conj(mapping_sequence(k+1))  -mapping_sequence(k+3);
                   mapping_sequence(k+3) -conj(mapping_sequence(k+1)) conj(mapping_sequence(k))  -mapping_sequence(k+2);
                   mapping_sequence(k+1) conj(mapping_sequence(k+3))  conj(mapping_sequence(k+2)) mapping_sequence(k) ;];
               X=X1*X2;  %酉矩陣調(diào)制    
                R=H*X;%+addtional_gaussian_noise(:,:,k);         %接收信號
                r1=-( conj(R(1,1))*R1(1,1)+R(1,2)*conj(R1(1,2))+R(1,3)*conj(R1(1,3))+conj(R(1,4))*R1(1,4) );   
                r2=-( conj(R(1,1))*R1(1,2)-R(1,2)*conj(R1(1,1))+R(1,3)*conj(R1(1,4))-conj(R(1,4))*R1(1,3) );   
                r3=-( conj(R(1,1))*R1(1,3)+R(1,2)*conj(R1(1,4))-R(1,3)*conj(R1(1,1))-conj(R(1,4))*R1(1,2) );  
                r4=-( conj(R(1,1))*R1(1,4)-R(1,2)*conj(R1(1,3))-R(1,3)*conj(R1(1,2))+conj(R(1,4))*R1(1,1) );  
                h=R1(1,1)*conj(R1(1,4))+conj(R1(1,1))*R1(1,4)-conj(R1(1,2))*R1(1,3)-R1(1,2)*conj(R1(1,3));    
                d=abs(R1(1,1))^2+abs(R1(1,2))^2+abs(R1(1,3))^2+abs(R1(1,4))^2;
                dh = [(1+j)/2 0;j/sqrt(2) 0;(-1+j)/2 0;-1/sqrt(2) 0;
                      (-1-j)/2 0;-j/sqrt(2) 0; (-1+j)/2 0;1/sqrt(2) 0;
                       0 exp(3j*pi/8)/sqrt(2);0 exp(5j*pi/8)/sqrt(2);
                       0 exp(7j*pi/8)/sqrt(2);0 exp(9j*pi/8)/sqrt(2);
                      0 exp(11j*pi/8)/sqrt(2);0 exp(13j*pi/8)/sqrt(2);
                      0 exp(15j*pi/8)/sqrt(2);0 exp(j*pi/8)/sqrt(2);];
                  
                 dc=[0 0 0 0;0 0 0 1;0 0 1 0;0 0 1 1;0 1 0 0;0 1 0 1;0 1 1 0;0 1 1 1;
                     1 0 0 0;1 0 0 1;1 0 1 0;1 0 1 1;1 1 0 0;1 1 0 1;1 1 1 0;1 1 1 1;];
                  for L=1:16
                %     A=[dh(L,1) 0 0 dh(L,2);0 conj(dh(L,1))  -conj(dh(L,2)) 0;0  -conj(dh(L,2)) conj(dh(L,1)) 0;dh(L,2) 0 0 dh(L,1)];
                %      D1(L)=real(trace(R'*R1*A));
                    D1(L)=real(-r1*dh(L,1)-r4*dh(L,2));
               %    D1(L)=d*(abs(dh(L,1))^2+abs(dh(L,2))^2)+2*real(r1*dh(L,1)+r4*dh(L,2)+h*(dh(L,1)*conj(dh(L,2))));
              end
              [minnumber position1]=max(D1);%求出最小的距離
                x1=dh(position1,1);  %估計出x1的值
                x4=dh(position1,2);  %估計出x2的值
                x14=dc(position1,:);

                for L=1:16
              %       B=[0 -conj(dh(L,1)) -conj(dh(L,2)) 0;dh(L,1) 0 0 -dh(L,2);dh(L,2) 0 0 -dh(L,1);0 conj(dh(L,2)) conj(dh(L,1)) 0];
              %       D2(L)=real(trace(R'*R1*B));
                    D2(L)=real(-r2*dh(L,1)-r3*dh(L,2));
               %   D2(L)=d*(abs(dh(L,1))^2+abs(dh(L,2))^2)+2*real(r2*dh(L,1)+r3*dh(L,2)+h*(dh(L,1)*conj(dh(L,2))));
                end
               
                     [minnumber position2]=max(D2);%求出最小的距離
                x2=dh(position2,1);  %估計出x1的值
                x3=dh(position2,2);  %估計出x2的值
                x23=dc(position2,:);
                if(X2(1,1)~=x1| X2(4,1)~=x4)  %計算Gray  mapping時bit錯誤的個數(shù)
             tot_err_bit((snr-min_snr)/2+1)=tot_err_bit((snr-min_snr)/2+1)+sum(abs(x14-bit_sequence(2*k-1:2*k+2))); 
                  end
                  
                 if(X2(2,1)~=x2 | X2(3,1)~=x3)  %計算Gray mapping時bit錯誤的個數(shù)
              tot_err_bit((snr-min_snr)/2+1)=tot_err_bit((snr-min_snr)/2+1)+sum(abs(x23-bit_sequence(2*k+3:2*k+6))); 
                 end
                X1=X; %參考酉矩陣
                R1=R; %參考信號    
           end
        end
    end        
    %在不同信噪比條件下，求系統(tǒng)SBER
    ber_snr=tot_err_bit./(tot_frame_tx * symbol_per_frame*2);
    figure(1);
   semilogy(min_snr:2:max_snr,ber_snr,'-o');
    grid
    xlabel('SNR(dB)');
    ylabel('Bit Error Ratio');
    legend('BER-SNR');