初识MIMO-OFDM（三）：完整的MIMO-OFDM仿真

零.缘起

MIMO-OFDM的仿真是一个大的不能再大的话题了，写得好的仿真平台是可以赚大钱的，我后续也会更新一个著名仿真平台的学习博客，不过我们这里先来学习一下陈老湿在通信家园看到的这个仿真https://zhuanlan.zhihu.com/p/392827532。这个博客主要是希望梳理每一个函数，然后给出2*2Alamouti+OFDM完整的仿真流程。
代码地址：https://github.com/liu-zongxi/MIMO_OFDM_simulation

发射机

一. 发射机步骤一：生成所需要发送的数据

%-----------------------生成比特-----------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年5月5日14点33分-----------------%
function Frame_bit = FrameBitGenerator(N_data, N_user, N_mod, N_symbol)
% 输入
% N_data:一个符号中要发送的符号数
% N_user:用户数
% N_mod:调制数
% N_symbol:一帧中有多少个符号
% 输出
% Frame_bit:一个（1，N_user）的结构体，包含data和num两个field，分别表示每个用户的数据和数据长度
    N_bit_per_symbol = N_mod * N_data;
    % 可以选择用户比特分布为：
    % 1) 所有用户比特相同。 
    N_bit_per_user_per_symbol = repmat(N_bit_per_symbol/N_user,1,N_user);

    % 2) 用户按照比例 (1:N_user)/((1 + N_user)*N_user/2) 发送比特数据
    % 目的是为了仿真不同用户的信息比特长度不同的情况,也可以修改得到其他的用户数据比例
    % N_bit_per_user_per_symbol = round([1:N_user]/((1 + N_user)*N_user/2) * N_bit_per_symbol);
    N_bit_per_user = N_bit_per_user_per_symbol * N_symbol;
    for iuser = 1:N_user
        % frame_bit_temp{iuser} = rand(N_bit_per_user(iuser), 1) > 0.5;
        Frame_bit(iuser).data = rand(N_bit_per_user(iuser), 1) > 0.5;
        Frame_bit(iuser).num = N_bit_per_user(iuser);
    end
end

该函数主要是生成所需要发送的数据，这里是为了扩展性用了结构体，实际上该仿真根本无法处理不同长度的数据。

一个OFDM符号所承载的比特数是他用于data传输的子载波数*调制数

由于四个用户等分数据，因此每个用户收到数据个数都要被除以N_user

最后我们得到的数据形式是

（N_data*N_mod/N_user*N_symbol, N_user）

二.发射机步骤二：子载波分配

%-----------------------子载波分配-----------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年5月5日16点35分-----------------%
function [index_data_per_user, Frame_zero_padding]=SubcarrierAllocation(Frame_bit_coded, index_data, N_data, N_user,N_symbol, N_mod, type_alloc)
Frame_zero_padding = cell(1,N_user);
index_data_per_user = cell(1,N_user);
% 计算填零的个数
N_subcarrier_per_user = N_data / N_user;    % 平均分配,每个用户所被分配的子载波个数
for iuser = 1:N_user
    N_zero_padding = N_mod * N_subcarrier_per_user * N_symbol - Frame_bit_coded(iuser).num;
    % 子载波分配方式，这里只给出相邻分配
    if type_alloc == "neighbour"
        index_data_per_user{iuser} = index_data((iuser-1)*N_subcarrier_per_user + 1: iuser*N_subcarrier_per_user)';
    end
    Frame_zero_padding{iuser} = [Frame_bit_coded(iuser).data; zeros(N_zero_padding, 1)];
end

该函数还是装神弄鬼了，他输出的是针对步骤一所得到的（N_data*N_mod/N_user*N_symbol, N_user）的数据，每一列的每一个symbol使用编号多少的子载波进行发送？

原仿真里给出了多种分配方式，但我这里单纯是实现一下，就是按照顺序分配的

最终输出的index_data_per_user是一个（N_data/N_user，N_user）的矩阵

Frame_zero_padding是一个零填充帧，如果现实情况中不是像我们步骤一中恰好生成那么多比特呢？这时要用0把他填满

三.发射机步骤三：调制

%-----------------------调制函数---------------------%
%-----------------------author:lzx--------------------------%
%-----------------------date:2022年5月5日23点22分-----------------%
function Frame_mod = Modulator(Frame_zero_padding, index_data_per_user, N_data, N_user, N_symbol, N_mod, N_subcarrier)

Frame_mod = zeros(N_subcarrier, N_symbol);
N_subcarrier_per_user = N_data / N_user;
for iuser = 1:N_user
    for isymbol = 1:N_symbol
        L_symbol = N_subcarrier_per_user* N_mod;
        index_symbol = (isymbol-1)*L_symbol+1 : isymbol*L_symbol;
        Symbol = Frame_zero_padding{iuser}(index_symbol);
        Symbol_premod = reshape(Symbol, N_mod, N_subcarrier_per_user);
        Symbol_mod = SymbolModulator(Symbol_premod);
        Frame_mod(index_data_per_user{iuser}, isymbol) = Symbol_mod.';
    end
end

%-----------------------MIMO-OFDM主函数---------------------%
%-----------------------author:lzx--------------------------%
%-----------------------date:2022年5月5日10点29分-----------------%
function Symbol_mod = SymbolModulator(Symbol_premod)

N_mod = size(Symbol_premod,1);

switch N_mod
    % BPSK调制
    case 1
        matrix_mapping = [-1 1];
        Symbol_mod = matrix_mapping(Symbol_premod + 1);
    % QPSK调制
    case 2
        %　比特的映射关系,00:-3/4*pi,01:3/4*pi,10: -1/4*pi,11: 1/4*pi
        matrix_mapping = exp(1j*[-3/4*pi 3/4*pi -1/4*pi 1/4*pi]);
        index_mod = [2 1]*Symbol_premod;
        %　把输入比特映射为符号
        Symbol_mod = matrix_mapping(index_mod + 1);
    % 8PSK调制
    case 3
        matrix_mapping = exp(1j*[0  1/4*pi 3/4*pi 1/2*pi  -1/4*pi -1/2*pi pi -3/4*pi ]);
        index_mod = [4 2 1]*Symbol_premod ;
        %　把输入比特映射为符号
        Symbol_mod = matrix_mapping(index_mod + 1);
    % 16QAM调制
    case 4
        % 映射关系参见说明文档
        m=1;
        for k=-3:2:3
            for l=-3:2:3
                % 对符号能量进行归一化
                matrix_mapping(m) = (k+1j*l)/sqrt(10);
            m=m+1;
            end
        end
        matrix_mapping = matrix_mapping([0 1 3 2 4 5 7 6 12 13 15 14 8 9 11 10]+1);
        index_mod = [8 4 2 1]*Symbol_premod ;
        Symbol_mod = matrix_mapping(index_mod + 1);
    % 64QAM调制         
    case    6
        % 映射关系参见说明文档
        m=1;
        for k=-7:2:7
            for l=-7:2:7
                % 对符号能量进行归一化
                matrix_mapping(m) = (k+1j*l)/sqrt(42); 
            m=m+1;
            end
        end
        matrix_mapping = matrix_mapping(...
         [[ 0  1  3  2  7  6  4  5]...
         8+[ 0  1  3  2  7  6  4  5]... 
         24+[ 0  1  3  2  7  6  4  5]...
         16+[ 0  1  3  2  7  6  4  5]...
         56+[ 0  1  3  2  7  6  4  5]...
         48+[ 0  1  3  2  7  6  4  5]...
         32+[ 0  1  3  2  7  6  4  5]...
         40+[ 0  1  3  2  7  6  4  5]]+1);
        index_mod = [32 16 8 4 2 1]*Symbol_premod ;
        Symbol_mod = matrix_mapping(index_mod + 1);
end

该函数是一个调制函数

分为两部分，首先在Frame选出一个user一个symbol所需要发送的比特，改变形状后进行调制，调制后放入步骤二所获得的对应的子载波中。

进行调制输入比特为（N_mod, N_subcarrier_per_user），输出一个（1，N_subcarrier_per_user），然后反转后放入到对应子载波中

调制采用的是非常巧妙的方式，首先获得格雷码的映射，然后一一对应，可惜仔细学习一下该函数。

四.步骤四：Alamouti编码

%-----------------------STBC时空块码-----------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年5月5日15点52分-----------------%
function Frame_STBC = STBCCoding(Frame_mod, N_subcarrier, N_symbol, N_Tx)
Frame_STBC = zeros(N_subcarrier, N_symbol, N_Tx);
if (mod(N_symbol,N_Tx))
    error('空时编码器输入符号不匹配,子程序st_coding出错');
else
    for ispace = 1:N_symbol/N_Tx
        if N_Tx == 2
            X1=Frame_mod(:,(ispace-1)*N_Tx+1);%取第一列，即第一个OFDM符号
            X2=Frame_mod(:,(ispace-1)*N_Tx+2);%取第二列，即第二个OFDM符号
            Symbol_STBC = [X1 X2;-conj(X2) conj(X1)];%alamouti编码
        elseif N_Tx == 4
            X1=Frame_mod(:,(ispace-1)*N_Tx+1);%取第一列，即第一个OFDM符号
            X2=Frame_mod(:,(ispace-1)*N_Tx+2);%取第二列，即第二个OFDM符号
            X3=Frame_mod(:,(ispace-1)*N_Tx+3);%取第一列，即第一个OFDM符号
            X4=Frame_mod(:,(ispace-1)*N_Tx+4);%取第二列，即第二个OFDM符号
            Symbol_STBC = [ X1 X2 X3  X4;...
                            -X2 X1 -X4 X3;...
                            -X3 X4 X1 -X2;...
                            -X4 -X3 X2 X1;...
                            conj(X1) conj(X2) conj(X3) conj(X4);...
                            -conj(X2) conj(X1) -conj(X4) conj(X3);...
                            -conj(X3) conj(X4) conj(X1) -conj(X2);...
                            -conj(X4) -conj(X3) conj(X2) conj(X1)];
        end
        for iant = 1:N_Tx
            Symbol_STBC_per_ant = reshape(Symbol_STBC(:,iant), N_subcarrier, N_Tx);
            Frame_STBC(:, (ispace-1)*N_Tx+1:ispace*N_Tx, iant) = Symbol_STBC_per_ant;
        end
    end
end

这是一个标准的STBC编码，每次取出相邻的两个symbol的符号（同时包含导频和数据，就是一个完整的symbol）（因为这才满足Alamouti编码0时刻和T时刻的要求）。然后进行STBC的编码把原本（N_subcarrier，2）变为（2*N_subcarrier，2）。之后对于编码后的每一列，是由一根天线发送的，把每一列reshape成（N_subcarrier， 2）（这个2是因为2*2编码），放入总矩阵中，当每一列都摆放后，得到一个（N_subcarrier，2，2）的矩阵，这其中第二个2是由于这个2是因为2*2编码，第三个2是由于两根天线，其实是一样的。

当所有的symbol都被取出后，最后得到一个（N_subcarrier，N_symbol，2）的矩阵，这就是编码结果了。

从结果上来看，它将原本（N_subcarrier，N_symbol）个信号生成为了（N_subcarrier，N_symbol，2）个信号，给两个天线分别发送，这就是发射分集。

五.步骤五：OFDM调制

%-----------------------调制函数---------------------%
%-----------------------author:lzx--------------------------%
%-----------------------date:2022年5月6日21点11分-----------------%
function Frame_transmit = OFDMModulator(Frame_pilot, N_sym, N_subcarrier, N_symbol, N_Tx, N_GI)

Frame_transmit = zeros(1,N_sym*N_symbol,N_Tx);

for iant = 1:N_Tx
    % ifft乘sqrt(N_subc)以保证变换前后能量不变
    % 我们假设频域的样点是在[-fs/2  fs/2]中的, fs是采样频率
    % 使用fftshift函数目的是使得变换前的频域样点转换到[0 fs]中,以满足IFFT变换的要求
    Frame_TD = sqrt(N_subcarrier) * ifft( fftshift( Frame_pilot(:,:,iant), 1 ) );
    Frame_CP = Frame_TD(N_subcarrier - N_GI + 1:N_subcarrier ,:);
    Frame_withCP = [Frame_CP; Frame_TD];
    % 转换为串行信号
    Frame_transmit(:,:,iant) = reshape( Frame_withCP, 1, N_sym*N_symbol);
end

该函数是对信号进行OFDM调制，调制按照天线来进行，将（N_subcarrier，N_symbol）大小的数据在列这一维上做fftshift然后进行ifft，这样还会得到（N_subcarrier，N_symbol）的数据，然后加CP，CP的大小是（N_GI, N_symbol），叠在数据的头上，最后，发送时转换为串行信号，发送一个（1，(N_GI+N_subcarrier)*N_symbol）的信号，多跟天仙的结果是（1，(N_GI+N_subcarrier)*N_symbol，2）.

这里为什么要做并串转换呢？实际上就是按顺序发送symbol，因为天线一次只能发送一个symbol呀~

信道

一.信道步骤一：加噪

这已经是老生常谈了

 Power_transmit = var(Frame_transmit);%发送信号功率
 N_noise = size(Frame_transmit,2);
 noise = NoiseGenerator(EbN0, Power_transmit, N_noise);
 Frame_noise = Frame_transmit+noise;

%-----------------------生成高斯噪声---------------------%
%-----------------------author:lzx--------------------------%
%-----------------------date:2022年5月6日22点43分-----------------%
function noise = NoiseGenerator(EbN0, Power_transmit, N_noise)
sigma = sqrt(Power_transmit/(2*EbN0));%标准差sigma
noise_real = randn(1,N_noise,2);
noise_real(:,:,1)=noise_real(:,:,1).*sigma(1);
noise_real(:,:,2)=noise_real(:,:,2).*sigma(2);
noise_imag = randn(1,N_noise,2);
noise_imag(:,:,1)=noise_imag(:,:,1).*sigma(1);
noise_imag(:,:,2)=noise_imag(:,:,2).*sigma(2);
noise=complex(noise_real,noise_imag);%复噪声序列

其实这个代码是比较愚蠢的，生成和发送信号等大的一个噪声矩阵，不详细说了

二.信道步骤二：信道

没有仿真的坐而论道，希望各位大佬指教

这个仿真里。。。。。居然没有信道模型，这其实是我最想知道的一个模型

目前对于信道的仿真，我们是割裂的。

在SISO-OFDM中，我们采用的是一个标准的TDL模型，他的大小是（1，Tau_max），然后进行卷积

在MIMO中，我们使用的是一个平坦信道，他是一个（N_Rx,NTx）的信道，直接进行相乘

那么在MIMO-OFDM中如何进行结合？因为每一个子载波都应该对应着一个子信道，信道的数量应该是(N_subcarrier, N_Tx*N_Rx);但这应该如何仿真呢？

接收机

首先要明确，接收机有N_user个

一. 步骤一：OFDM解调

%-----------------------OFDM的接收---------------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年5月5日16点35分-----------------%
function Frame_recieve = OFDMDemodulator(Frame_noise, N_sym, N_subcarrier,N_symbol,N_Rx, N_GI)
Frame_recieve = zeros(N_subcarrier,N_symbol,N_Rx);

for iant = 1:N_Rx
    Frame_symbol = reshape(Frame_noise(1,:,iant), N_sym, N_symbol);
    Frame_noGI = Frame_symbol(N_GI+1:end, :);
    % fft乘1/sqrt(N_subc)以保证变换前后能量不变
    % 我们假设频域的样点是在[-fs/2  fs/2]中的, fs是采样频率
    % fftshift目的是使得变换后的频域样点在[-fs/2  fs/2]中,而不是[0 fs]中
    Frame_recieve(:,:,iant) = fftshift(1/sqrt(N_subcarrier) * fft( Frame_noGI ), 1);
end
% Frame_recieve = Frame_FD(:, 1:N_symbol ,:);     % 数据OFDM符号, 包括导频符号

作为天线，自然应该按天线处理

接收到的信号重新串并转换为（N_sym, N_symbol)，因为我们肯定是一个symbol一个symbol接收到的，去掉CP然后反做fft即可

最后得到的输出是（N_subcarrier，N_symbol），按天线接收后，得到的是（N_subcarrier，N_symbol，2）不过我这里存在疑问，接收端如何能分开发射端两根天线？我对此感到不解

二. 接收机二：STBC解码

%-----------------------STBC解码---------------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年5月7日10点07分-----------------%
function Frame_decoded = STBCDecoding(Frame_recieve, H, N_subcarrier, N_Tx,N_Rx, N_symbol)
Frame_decoded = zeros(N_subcarrier, N_symbol);
% 把H_freq转化为空时译码器的输入格式,为一个N_subc*N_ant_pair的矩阵,每列表示:
% 1-->1 ,1-->2,...,1-->N_Rx_ant, ... ,N_Tx_ant-->1, N_Tx_ant-->2,..., N_Tx_ant-->N_Rx_ant
if (N_Tx == 2)&&(N_Rx == 2)
    for ispace = 1:N_symbol/N_Tx
        % 构造输入进空时译码器的符号,用t表示时间号,a表示天线号,其格式为:
        % [Recv(t1,a1) Recv(t2,a1) Recv(t1,a2) Recv(t2,a2)].
        R = [];
        for iant = 1:N_Rx
            R = [R  Frame_recieve(:,(ispace-1)*N_Tx+1:ispace*N_Tx,iant) ];%对于2X2MIMO来说，R矩阵即为[Xe  -conj(Xo)  Xo  conj(Xe)]
        end
        H11=H(:,1);%1发送天线--》1接收天线信道参数
        H12=H(:,2);%1发送天线--》2接收天线信道参数
        H21=H(:,3);%2发送天线--》1接收天线信道参数
        H22=H(:,4);%2发送天线--》2接收天线信道参数
        R11=R(:,1);%1发送天线发送的第一个符号
        R12=R(:,2);%1发送天线发送的第二个符号
        R21=R(:,3);%2发送天线发送的第一个符号
        R22=R(:,4);%2发送天线发送的第二个符号
        for i=1:1:N_subcarrier
            X1(i,1)=(R11(i)*conj(H11(i))+conj(R12(i))*H21(i)+R21(i)*conj(H12(i))+conj(R22(i))*H22(i))/( H11(i)*conj(H11(i)) + H21(i)*conj(H21(i)) + H12(i)*conj(H12(i)) + H22(i)*conj(H22(i)));
            X2(i,1)=(R11(i)*conj(H21(i))-conj(R12(i))*H11(i)+R21(i)*conj(H22(i))-conj(R22(i))*H12(i))/( H11(i)*conj(H11(i)) + H21(i)*conj(H21(i)) + H12(i)*conj(H12(i)) + H22(i)*conj(H22(i)));
        end
        Sybmol_decoded = 2*[X1 X2]; % 非常扯淡，因为上面的分母是错误的
        Frame_decoded(:,(ispace-1)*N_Tx+1:ispace*N_Tx) = Sybmol_decoded;
    end
end

首先我要申明，我认为该仿真是错的Sybmol_decoded = 2*[X1 X2]; % 非常扯淡，因为上面的分母是错误的为什么莫名其妙的×2呢，因为上面( H11(i)*conj(H11(i)) + H21(i)*conj(H21(i)) + H12(i)*conj(H12(i)) + H22(i)*conj(H22(i)));是错误的，应该只选取其中两个。代码的实现更是丑陋，所以我也没有重写的欲望了。她也没有解决我想知道的信道的问题，

对于每两列进行一个操作，再加上两根天线，生成会原本的（N_subcarrier*2,2）的矩阵，然后解算。的总之他会把原本（N_subcarrier，N_symbol，2）的信号接收后得到N_subcarrier，N_symbol).

三.接收机步骤三：解调

%-----------------------QAM解调-----------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年5月7日11点51分-----------------%
function Frame_demod = Demodulator(Frame_decoded, index_data_iuser, N_mod, N_symbol)

Frame_demod = [];

for isymbol = 1:N_symbol
        Symbol_demod = SymbolDemodulator(Frame_decoded(index_data_iuser, isymbol).', N_mod);
        Frame_demod = [Frame_demod; Symbol_demod(:)];
end

%-----------------------每个符号的QAM解调-----------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年5月7日11点51分-----------------%
function Symbol_demod = SymbolDemodulator(Symbol, N_mod)

switch N_mod

    % BPSK解调
    case 1
        Symbol_demod = real(Symbol) > 0; 

    % QPSK解调
    case 2 
        % 由QPSK的星座图可以观察到
        bit0 = real(Symbol) ; 
        bit1 = imag(Symbol) ;
        
        % 得到2行, 列数为符号数的输出矩阵
        Symbol_demod(1,:) = bit0 > 0;
        Symbol_demod(2,:) = bit1 > 0;

    % 8PSK解调
    case 3
        % 参见8PSK的星座图
        bit0 = -imag( Symbol * exp(1j*pi/8)) ;
        % bit1和bit2解调,都需要进行星座旋转
        bit1 = -real(Symbol * exp(1j*pi/8)) ;
        
        bit2 = [];
        for isymbol = 1:length(Symbol)
            tmp = Symbol(isymbol) * exp(-1j*pi/8); 
            if ((real(tmp) <0) && (imag(tmp) >0)) || ((real(tmp) >0) && (imag(tmp) <0))
                bit2 = [bit2 0];
            else
                bit2 = [bit2 1];
            end   
        end
        
        Symbol_demod(1,:) = bit0 >0;
        Symbol_demod(2,:) = bit1 >0;
        Symbol_demod(3,:) = bit2 ;    %　已经硬判决
            % 16QAM解调   
    case    4
        bit0 = real(Symbol);
        bit2 = imag(Symbol);

        % 以bit1的生成来说明方法:
        % 2/sqrt(10) 为临界值, abs(real(sym))大于此, 则bit1为负,硬判决得到0 ; 反之为正
        bit1 = 2/sqrt(10)-(abs(real(Symbol)));
        bit3 = 2/sqrt(10)-(abs(imag(Symbol)));

        Symbol_demod(1,:) = bit0 > 0;
        Symbol_demod(2,:) = bit1 > 0;
        Symbol_demod(3,:) = bit2 > 0;
        Symbol_demod(4,:) = bit3 > 0;
            % 64QAM解调         
    case    6       
        bit0 = real(Symbol);
        bit3 = imag(Symbol);
        bit1 = 4/sqrt(42)-abs(real(Symbol));
        bit4 = 4/sqrt(42)-abs(imag(Symbol));
        for m=1:size(Symbol,2)
            for k=1:size(Symbol,1)
                if abs(4/sqrt(42)-abs(real(Symbol(k,m)))) <= 2/sqrt(42)  
                    bit2(k,m) = 2/sqrt(42) - abs(4/sqrt(42)-abs(real(Symbol(k,m))));
                elseif abs(real(Symbol(k,m))) <= 2/sqrt(42) 
                    bit2(k,m) = -2/sqrt(42) + abs(real(Symbol(k,m)));
                else
                    bit2(k,m) = 6/sqrt(42)-abs(real(Symbol(k,m)));
                end
      
                if abs(4/sqrt(42)-abs(imag(Symbol(k,m)))) <= 2/sqrt(42)  
                    bit5(k,m) = 2/sqrt(42) - abs(4/sqrt(42)-abs(imag(Symbol(k,m))));
                elseif abs(imag(Symbol(k,m))) <= 2/sqrt(42) 
                    bit5(k,m) = -2/sqrt(42) + abs(imag(Symbol(k,m)));
                else
                    bit5(k,m) = 6/sqrt(42)-abs(imag(Symbol(k,m)));
                end
            end
        end

        Symbol_demod(1,:) = bit0 > 0;
        Symbol_demod(2,:) = bit1 > 0;
        Symbol_demod(3,:) = bit2 > 0;
        Symbol_demod(4,:) = bit3 > 0;
        Symbol_demod(5,:) = bit4 > 0;
        Symbol_demod(6,:) = bit5 > 0;
end

这里对于QAM，我没有自己写过就不班门弄斧了，但我大致知道他的思路

解调时，也是对于每一个symbol的每一个用户的那一部分子载波进行解调，挑选出属于该接收机 子载波然后全部并联起来，最后得到一个（N_data*Nysmbol, 1）便于计算误码率。

总结

总的来说，这个仿真我觉得不行!而且我认为OFDM的发送解调和STBC的发送解调顺序可能有问题。如果后续看到更好的仿真，我会继续更新的！