2mimo.m

在MIMO系统中的二乘二系统的模型代码

%**************** 内容：发射端已知SVD ****************
%****************       4*4 STBC U 1  ****************
%**************** 编程人：徐娜 ***********************
%**************** 时间：2008.2.27 ********************

clear;
clc;

%**************************** 初始化 ******************************

Num_Tr = 4;                                      % 发射天线数
Num_Re = 4;                                      % 接收天线数
Len_Frame = 100;                                 % 每帧数据长度
Data = zeros(1,Num_Tr*Len_Frame);                % 原始数据(分层空时码原始数据是矩阵，这里是串行数据向量，数据个数与分层一样)
Data_Modulated=zeros(1,Num_Tr*Len_Frame);        % 已调信号
Transmit = zeros(Num_Tr,Num_Tr*Len_Frame*2);     % 发射符号矩阵
Receive = zeros(Num_Re,Num_Tr*Len_Frame*2);      % 加噪声时的接收信号

SNR_Begin = 0;                                   % SNR单位为dB
SNR_Step = 5;
SNR_End = 30;
Step_Num = (SNR_End - SNR_Begin)/SNR_Step+1;     % SNR坐标的个数
Err = zeros(Step_Num,1);                         % 对各SNR下的误码个数初始化，以便在Monte Carlo仿真时叠加

a = 0;
for SNR0 = SNR_Begin : SNR_Step : SNR_End        % 将信噪比的取值点存成向量 
    a = a+1; 
    SNR(a,1) = SNR0;
end

table=exp(j*[pi/4 3*pi/4 -3*pi/4 -pi/4]);                      % 产生QPSK信号星座图模板                
        temp=[];
        temp1=[];
        temp2=[];
        for a=1:4
            for b=1:4
                for c=1:4
                    for d=1:4
                        temp1=[table(a);table(b);table(c);table(d)];
                        temp2=[temp2 temp1];
                        temp1=[];
                    end
                    temp=[temp temp2];
                    temp2=[];
                end
                temp=[temp temp2];
                temp2=[];
            end
            temp=[temp temp2];
            temp2=[];
        end                                                            % 模板矩阵temp的维数为4*256 

t = waitbar(0,'Please wait...');

Num_Mon = 2500;                               % Monte Carlo仿真次数(发送帧数)

%************************* Monte Carlo仿真 *************************

for MM = 1:Num_Mon

%**************************** 发射信号 *****************************

Fd = 1;                             % 消息序列的采样速率
Fs = 1*Fd;                          % 已调信号的采样速率
Md = 4;                             % 调制进制数,如: 若Md=2,则Md_PSK=BPSK;若Md=4,则Md_PSK=QPSK

Data = randint(1,Num_Tr*Len_Frame,Md);               % 随机产生原始信号(4QAM调制，则此处Data为0-3的整数)
Data_Modulated = dmodce(Data,Fd,Fs,'QAM',Md);        % 生成Md进制的QAM调制信号

for k=1:Len_Frame;                                   %产生发射信号
    Transmit(1,8*k-7)= Data_Modulated(1,4*k-3);
    Transmit(1,8*k-6)=-Data_Modulated(1,4*k-2);
    Transmit(1,8*k-5)=-Data_Modulated(1,4*k-1);
    Transmit(1,8*k-4)=-Data_Modulated(1,4*k);
    Transmit(1,8*k-3)= conj(Data_Modulated(1,4*k-3));
    Transmit(1,8*k-2)=-conj(Data_Modulated(1,4*k-2));
    Transmit(1,8*k-1)=-conj(Data_Modulated(1,4*k-1));
    Transmit(1,8*k)  =-conj(Data_Modulated(1,4*k));
    Transmit(2,8*k-7)= Data_Modulated(1,4*k-2);
    Transmit(2,8*k-6)= Data_Modulated(1,4*k-3);
    Transmit(2,8*k-5)= Data_Modulated(1,4*k);
    Transmit(2,8*k-4)=-Data_Modulated(1,4*k-1);
    Transmit(2,8*k-3)= conj(Data_Modulated(1,4*k-2));
    Transmit(2,8*k-2)= conj(Data_Modulated(1,4*k-3));
    Transmit(2,8*k-1)= conj(Data_Modulated(1,4*k));
    Transmit(2,8*k)  =-conj(Data_Modulated(1,4*k-1));
    Transmit(3,8*k-7)= Data_Modulated(1,4*k-1);
    Transmit(3,8*k-6)=-Data_Modulated(1,4*k);
    Transmit(3,8*k-5)= Data_Modulated(1,4*k-3);
    Transmit(3,8*k-4)= Data_Modulated(1,4*k-2);
    Transmit(3,8*k-3)= conj(Data_Modulated(1,4*k-1));
    Transmit(3,8*k-2)=-conj(Data_Modulated(1,4*k));
    Transmit(3,8*k-1)= conj(Data_Modulated(1,4*k-3));
    Transmit(3,8*k)  = conj(Data_Modulated(1,4*k-2));
    Transmit(4,8*k-7)= Data_Modulated(1,4*k);
    Transmit(4,8*k-6)= Data_Modulated(1,4*k-1);
    Transmit(4,8*k-5)=-Data_Modulated(1,4*k-2);
    Transmit(4,8*k-4)= Data_Modulated(1,4*k-3);
    Transmit(4,8*k-3)= conj(Data_Modulated(1,4*k));
    Transmit(4,8*k-2)= conj(Data_Modulated(1,4*k-1));
    Transmit(4,8*k-1)=-conj(Data_Modulated(1,4*k-2));
    Transmit(4,8*k)  = conj(Data_Modulated(1,4*k-3));
end

if Md == 2                                         % 每根发送天线发送功率归一化
   Transmit = Transmit;                            % 此处只适用Md=2和4，对于16QAM等，还需增加代码 
elseif Md == 4
   Transmit = Transmit/sqrt(2);
else 
   Md
   disp('没有定义此时的情况，请修改程序');
return
end

%**************************** MIMO信道 *****************************

%H = zeros(Num_Re,Num_Tr);                     %Rayleigh信道
%x = randn(Num_Re,Num_Tr);
%y = randn(Num_Re,Num_Tr);
%H(:,:) =  x + j*y ;

%H =random('Rayleigh',sqrt(0.5),Num_Re,Num_Tr)+i*random('Rayleigh',sqrt(0.5),Num_Re,Num_Tr);     %Rayleigh信道 

H =random('Rayleigh',sqrt(0.5),Num_Re,Num_Tr)+i*random('Rayleigh',sqrt(0.5),Num_Re,Num_Tr);
[Q R]=qr(H);                          %得到酉矩阵
H_U=Q;                              %酉矩阵的奇异值全为1，乘以4将奇异值改为4
 
%************************ 接收信号(无噪声) *************************

Receive_noN = zeros( Num_Re, Num_Tr*Len_Frame*2 );     % 不加噪声时的接收信号
Receive_noN = H_U * Transmit;                            % 接收信号

%********************* 不同SNR条件下接收译码 ***********************

a = 0;
for SNR0 = SNR_Begin : SNR_Step : SNR_End            % 将信噪比的取值点存成向量 
    a = a+1; 
    
    % ---------------------- 定义噪声 -------------------------
    
    Power_Sig = sum( sum( abs(H_U).^2 ) );             % 所有接收天线上总的信号功率
    Power_Noise = Power_Sig * 10^(-SNR0/10);         % 采用等增益合并方式，各接收天线上噪声功率相等
                                                     % 通过改变(降低)每次的噪声功率来改变(增加)信噪比   
    Noise = 1/sqrt(2) * ( normrnd( 0,sqrt(Power_Noise),Num_Tr,(Num_Tr*Len_Frame*2) ) + j*normrnd( 0,sqrt(Power_Noise),Num_Tr,(Num_Tr*Len_Frame*2) )  ); % 定义噪声 
                                                     % 除以sqrt(2)是因为复高斯噪声功率为实部虚部平方和
    
    % ---------------------- 接收信号 -------------------------
    
    Receive = Receive_noN + Noise;

    % ------------------------ 译码 ---------------------------
    
    Receive_Decode=zeros(Num_Re,Len_Frame);         % 书94页
    for k=1:Len_Frame
        for p=1:Num_Re
            Receive_Decode(1,k)=Receive_Decode(1,k)+Receive(p,8*k-7)*conj(H(p,1))+Receive(p,8*k-6)*conj(H(p,2))+Receive(p,8*k-5)*conj(H(p,3))+Receive(p,8*k-4)*conj(H(p,4))+conj(Receive(p,8*k-3))*H(p,1)+conj(Receive(p,8*k-2))*H(p,2)+conj(Receive(p,8*k-1))*H(p,3)+conj(Receive(p,8*k))*H(p,4);
            Receive_Decode(2,k)=Receive_Decode(2,k)+Receive(p,8*k-7)*conj(H(p,2))-Receive(p,8*k-6)*conj(H(p,1))-Receive(p,8*k-5)*conj(H(p,4))+Receive(p,8*k-4)*conj(H(p,3))+conj(Receive(p,8*k-3))*H(p,2)-conj(Receive(p,8*k-2))*H(p,1)-conj(Receive(p,8*k-1))*H(p,4)+conj(Receive(p,8*k))*H(p,3);
            Receive_Decode(3,k)=Receive_Decode(3,k)+Receive(p,8*k-7)*conj(H(p,3))+Receive(p,8*k-6)*conj(H(p,4))-Receive(p,8*k-5)*conj(H(p,1))-Receive(p,8*k-4)*conj(H(p,2))+conj(Receive(p,8*k-3))*H(p,3)+conj(Receive(p,8*k-2))*H(p,4)-conj(Receive(p,8*k-1))*H(p,1)-conj(Receive(p,8*k))*H(p,2);
            Receive_Decode(4,k)=Receive_Decode(4,k)+Receive(p,8*k-7)*conj(H(p,4))-Receive(p,8*k-6)*conj(H(p,3))+Receive(p,8*k-5)*conj(H(p,2))-Receive(p,8*k-4)*conj(H(p,1))+conj(Receive(p,8*k-3))*H(p,4)-conj(Receive(p,8*k-2))*H(p,3)+conj(Receive(p,8*k-1))*H(p,2)-conj(Receive(p,8*k))*H(p,1);
        end
    end
    h=norm(H_U,'fro');
    Receive_Decode=Receive_Decode/(2*(h.^2));
    
        %----------------- STBC的ML译码 -------------------               
                                              
        dif=[];
        Decode=[];                                                     % 译码后的结果
        for m=1:Len_Frame
            for n=1:256
                dif(:,n)=abs(Receive_Decode(:,m)-temp(:,n)).^2;
            end
            f=sum(dif);
            [y1 x1]=min(f);
            Decode(:,m)=[temp(:,x1)];
            dif=[];
            f=[];
        end                                                          
        
        
        
        % -------------- 求误码个数(不是误bit个数) ----------------
        
        for q=1:Len_Frame
            Data_Matrix(1,q)=Data(1,4*q-3);
            Data_Matrix(2,q)=Data(1,4*q-2);
            Data_Matrix(3,q)=Data(1,4*q-1);
            Data_Matrix(4,q)=Data(1,4*q);
        end
            
        Data_Est = ddemodce(Decode,Fd,Fs,'QAM',Md);       % 对译码结果进行判决，并转化为 0-3
        Err(a,1) = nnz(Data_Matrix - Data_Est) + Err(a,1);

end	    % 结束不同SNR的循环

waitbar(MM/Num_Mon);

end     % 结束不同 Monte Carlo 循环

close(t);

SER=Err/(Num_Tr*Len_Frame*Num_Mon)

%************************** 显示误码率曲线 **************************

semilogy(SNR,SER,'o-');
xlabel('SNR(dB)');
ylabel('SER');