初识MIMO(五)：CSI反馈及其仿真

零代码地址

https://github.com/liu-zongxi/MIMO_simulation

一. 发射端的信道估计

这章其实给了我隔靴搔痒的感觉，CSI的估计要尽可能快的完成，这里给出了这么几条路

利用信道的互异性，这样可以在发射端估计
量化CSI ，就是把CSI压缩后传输
利用码本，码本表示CSI

具体怎么做呢，书中没有给出，以后慢慢读论文吧！

二. CSI 在接收端的利用技术——Alamouti预编码及其仿真

1. 概念

2.代码展示

%-----------------------Alamouti预编码对比-------------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年4月30日09点58分-----------------%
%% 设置参数
clear;clf;
L_frame = 400; % 帧长度和帧个数
N_iter = 1500;
Nmod = 2;
M = 2^Nmod;SNRs_dB = 0:2:20;       %信噪比
SNRs = 10.^(SNRs_dB / 10.);
N_SNR = length(SNRs_dB);T = 2;        % Alamouti长度T
N_M = 2;      % 列向量长度，码字长度
L_code = 64;% 码本长度
NT = 4;     % 发射和接收天线个数
NR = 1;
N_frame_bit = NT*Nmod*L_frame;  %一帧比特数
N_totoal_bit = N_frame_bit*N_iter;% 总比特数
code_book = CodebookGenerator(NT, N_M, L_code);%生成码本 y = zeros(2,1,L_frame);
Rx = zeros(1,2,L_frame);
BERs = zeros(2, N_SNR);fprintf('====================================================\n');
fprintf('  Precoding transmission');
fprintf('\n  %d x %d MIMO\n  %d QAM', NT,NR,M);
fprintf('\n  Simulation bits : %d', N_totoal_bit);
fprintf('\n====================================================\n');%% 主函数
for isnr = 1:N_SNRSNR = SNRs(isnr);sigma = sqrt(NT/(2*SNR));n_biterror = 0;for iiter = 1:N_iterframe_origin = randi([0,1], L_frame, Nmod*N_M); % 生成随机数% 调制frame_mod = QPSKMod(frame_origin,L_frame,N_M);% Alamouti发送frame_transmit = reshape(frame_mod, N_M, 1, L_frame);frame_alamouti = [frame_transmit(1,1, :) -conj(frame_transmit(2,1, :));frame_transmit(2,1, :) conj(frame_transmit(1,1, :))];% 信道和噪声Hiid = (randn(NR,NT)+1j*randn(NR,NT))/sqrt(2);arg_W = zeros(1, L_code);for ih = 1:L_codearg_W(ih) = norm(Hiid*code_book(:,:,ih),'fro');end[arg_val, arg_Index] = max(arg_W);HW = Hiid*code_book(:,:,arg_Index);for iframe = 1:L_frameRx(:,:,iframe) = HW*frame_alamouti(:,:,iframe)+sigma*(randn(NR,T)+1j*randn(NR,T))end% 接收HW_re = [conj(HW(1)) HW(2); conj(HW(2)) -HW(1)];Rx_re = [Rx(1,1,:); conj(Rx(1,2,:))];for iframe = 1:L_framey(:,:,iframe) = HW_re* Rx_re(:,:,iframe)./norm(HW)^2;endframe_re = reshape(y, L_frame, N_M);frame_demod = QPSKDemod(frame_re,L_frame,N_M);% 计算误码率n_biterror_tmp = sum(sum(abs(frame_demod - frame_origin)));n_biterror = n_biterror + n_biterror_tmp;endBERs(1, isnr) = n_biterror/(N_iter*L_frame*Nmod*N_M);
end
% 4*1 无编码
BERs(2,:) = [0.345555000000000  0.245437500000000   0.151757500000000   0.0813550000000000  0.0365075000000000  0.0133750000000000  0.00412000000000000 0.00103500000000000 0.000242500000000000    5.00000000000000e-05    2.25000000000000e-05];
semilogy(SNRs_dB,BERs(1,:),"-^", SNRs_dB,BERs(2,:),'-+');
axis([SNRs_dB([1 end]) 1e-6 1e0]);
xlabel('SNR[dB]'), ylabel('BER');
legend('Precoded Alamouti','No Precoded Alamouti');

%-----------------------码本生成-------------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年4月30日09点58分-----------------%
function [code_book]=CodebookGenerator(NT, M, L)
%输入
%NT: 天线数
%M: 列向量长度，看书码字长度
%L: 码本长度
%输出
% code_book：码本
if NT == 2 && M == 1 && L ==8cloumn_index = [1];rotation_vector = [1,0];
elseif NT == 3 && M == 1 && L ==32cloumn_index=[1]; rotation_vector=[1 26 28];
elseif NT == 4 && M == 2 && L ==32cloumn_index=[1 2]; rotation_vector=[1 26 28];
elseif NT == 4 && M ==1 && L == 64cloumn_index=[1]; rotation_vector=[1 8 61 45];
elseif NT == 4 && M ==2 && L == 64cloumn_index=[0 1]; rotation_vector=[1 7 52 56];
elseif NT == 4 && M ==3 && L == 64cloumn_index=[0 2 3]; rotation_vector=[1 8 61 45];
elseerror("没有对应的码本");
end
index_w = 0:NT-1;
w = exp(1j*2*pi/NT*(index_w).'*index_w)/sqrt(NT);
w_1 = w(:, cloumn_index+1);
theta = diag(exp(1j*2*pi/L*rotation_vector));
code_book(:,:,1) = w_1;
for i = 1:L-1code_book(:,:,i+1) = theta*code_book(:,:,i);
end

3.对于代码的一些思考

原书中对于代码的介绍简直就是胡扯蛋！

我们先来看一下原书对于此段代码的论述

这简直是一派胡言啊，你所用的码本和你仿真天线个数根本就无法对应啊，这如何能对应的上呢？

在我修改后我对代码是这样理解的。首先码本引入了一个非常重要的中间量M，在这段仿真中，M=2，这意味着我们使用的Alamouti编码应该是2天线的，不过码本W有能力把他转换为4天线发送，这样发送同时可以抵消H的影响，让性能更好。

所以我们最后比较也是使用了4*1的Alamouti进行比较的

在这次代码中，我的仿真终于是严格按照书上公式来的了

之前的代码由于历史遗留问题，一直是二维矩阵的形式，这次改为了循环+三维矩阵的形式，虽然运行时间长了，但是更容易理解了。

发送的x应该是一个 (NT, 1, L_frame)的三维矩阵

H则是一个（NR,NT）的二维矩阵

noise是一个（NR，T）的二维矩阵，这里的T代表的是Alamouti编码里的T

接收到的y是一个（NR,1,L_frame）的三维矩阵

噪信比

sigma = sqrt(NT/(2*SNR));

之前一直迷迷糊糊的仿真

天线越多，每根天线的功率越小，sigma也就越大

三. CSI 在接收端的利用技术——ZF及MMSE预编码及其仿真

1.概念

2.代码

%---------------------ZF 和 MMSE预编码----------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:17点18分-----------------%
%% 参数设置
clear; clf;
NT = 4;
NR = 4;     % 天线数
L_frame = 1000;  %帧长度
N_iter = 1000;  % 循环次数
SNRs_dB = 3:1:30;   % 信噪比
SNRs = 10.^(SNRs_dB./10);
N_SNR = length(SNRs);
Nmod = 2;       % QPSK
N_case = 4;          % 不同类型
BERs = zeros(N_case, N_SNR);
gss = ["-kx" "-^" "-ro" "-b>" "-g<" "-m+"];   % 画图图像，注意使用双引号
y = zeros(NR,1,L_frame);
x_tilde = zeros(NR,1,L_frame);
%% 主函数
for icase = [1 2 3 4]gs = gss(icase);if icase == 1% ZFW_pre_formula = @(Hiid, sigma, NT) eye(size(Hiid));W_formula = @(Hiid, sigma, NT) inv(Hiid'*Hiid) * Hiid';% W_formula = @(Hiid, sigma, NT) inv(Hiid);beta_formula = @(Hiid, sigma, NT) 1;elseif icase == 2% MMSEW_pre_formula = @(Hiid, sigma, NT) 1;W_formula = @(Hiid, sigma, NT) inv(Hiid'*Hiid +sigma.^2*diag(ones(1,NT))) * Hiid';beta_formula = @(Hiid, sigma, NT) 1;elseif icase == 3% W_pre_formula = @(Hiid, sigma, NT) sqrt(NT/trace(inv(Hiid)*inv(Hiid)'))*inv(Hiid);W_pre_formula = @(Hiid, sigma, NT) sqrt(NT/trace((Hiid' * inv(Hiid*Hiid'))*(Hiid' * inv(Hiid*Hiid'))'))*Hiid' * inv(Hiid*Hiid');W_formula = @(Hiid, sigma, NT) eye(size(Hiid));beta_formula = @(Hiid, sigma, NT) sqrt(NT/trace((Hiid' * inv(Hiid*Hiid'))*(Hiid' * inv(Hiid*Hiid'))'));elseif icase == 4W_pre_formula = @(Hiid, sigma, NT) sqrt(NT/trace((Hiid'*inv(Hiid*Hiid'+sigma.^2*eye(NR,NR)))*(Hiid'*inv(Hiid*Hiid'+sigma.^2*eye(NR,NR)))'))*Hiid'*inv(Hiid*Hiid'+sigma.^2*eye(NR,NR));W_formula = @(Hiid, sigma, NT) eye(size(Hiid));beta_formula = @(Hiid, sigma, NT) sqrt(NT/trace((Hiid'*inv(Hiid*Hiid'+sigma.^2*eye(NR,NR))*(Hiid'*inv(Hiid*Hiid'+sigma.^2*eye(NR,NR)))')));endfor isnr = 1:N_SNRSNR = SNRs(isnr);n_biterror = 0;for iiter = 1:N_iter% 生成数据frame_origin = randi([0,1],L_frame,Nmod*NT);% QPSK调制frame_mod=QPSKMod(frame_origin,L_frame, NT);% 调整为标准形式frame_transmit = reshape(frame_mod, NT, 1, L_frame);% 生成信道，SIMO有NR个信道Hiid = (randn(NR,NT)+1j*randn(NR,NT))./sqrt(2);%  AWGN噪声sigma = sqrt(NT/(2*SNR));noise = sigma*(randn(NR, 1) + 1j*randn(NR, 1));% 预编码矩阵W_pre = W_pre_formula(Hiid, sigma, NT);% 接收信号beta = beta_formula(Hiid, sigma, NT);for iframe = 1:L_frame
%                 if icase == 1
%                     y(:,:,iframe) = (Hiid*frame_transmit(:,:,iframe)+noise)/beta;
%                 endy(:,:,iframe) = (Hiid * W_pre * frame_transmit(:,:,iframe)+noise)/beta;end% 信号检测W = W_formula(Hiid, sigma, NT);for iframe = 1:L_framex_tilde(:,:,iframe) = W * y(:,:,iframe);end% 解调frame_re = reshape(x_tilde, L_frame, NT);frame_demod = QPSKDemod(frame_re,L_frame,NT);% 计算误码率n_biterror_tmp = sum(sum(abs(frame_demod - frame_origin)));n_biterror = n_biterror + n_biterror_tmp;endBERs(icase, isnr) = n_biterror/(N_iter*L_frame*Nmod*2);end
end
%% 画图
semilogy(SNRs_dB,BERs(1,:),gss(1), SNRs_dB,BERs(2,:),gss(2), SNRs_dB,BERs(3,:),gss(3), SNRs_dB,BERs(4,:),gss(4));
axis([SNRs_dB([1 end]) 1e-6 1e0])
xlabel('SNR[dB]'), ylabel('BER');
legend('ZF','MMSE', "Pre-ZF", "Pre-MMSE");

3.一些思考

预编码的MMSE和ZF想法是很简单的，均衡就是在接收端 $H^{-1}Hx$ ,而预编码就是 $HH^{-1}x$ ，由于不受到噪声影响，就会效果好

不过我有一个疑惑，应该beta大于1才能减小噪声的干扰吗？为什么都可以，希望大佬解答

四.天线选择技术

1.概念

2.代码

%-----------------------Alamouti预编码对比-------------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年4月30日09点58分-----------------%
%% 参数设置
clear;
clf;
NT=4;
NR=4;
N_iter=1000;
I=eye(NR,NR);
gss=['-kx';'-ro';'-k^';'-g+'];
SNRs_dB=0:2:20;
SNRs = 10.^(SNRs_dB/10);
N_SNR = length(SNRs);
Capacitys = zeros(1,N_SNR);
%% 主函数
for n_ant=1:NTif n_ant>NT || n_ant<1error('n_ant must be between 1 and NT!');elseTypes = nchoosek([1:NT],n_ant);N_type = size(Types,1);endfor isnr=1:N_SNRSNR = SNRs(isnr)/n_ant;
%       rand('seed',1); randn('seed',1);  sum_capacity = 0;for iiter=1:N_iterH = (randn(NR,NT)+1j*randn(NR,NT))/sqrt(2);log_SH = zeros(1,N_type);for itype=1:N_typeH_sel = H(:,Types(itype,:)); log_SH(itype)=log2(real(det(I+SNR*H_sel*H_sel'))); % Eq.(12.22)endsum_capacity = sum_capacity + max(log_SH);endCapacitys(isnr) = sum_capacity/N_iter;endplot(SNRs_dB,Capacitys,gss(n_ant,:), 'LineWidth',2); hold on;
end
xlabel('SNR[dB]'), ylabel('bps/Hz'), grid on;
legend('sel-ant=1','sel-ant=2','sel-ant=3','sel-ant=4')

%-----------------------降低复杂度的天线选择-------------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年4月30日09点58分-----------------%
%% 参数设置
clear; clf;
method_sel = 1; % 0/1 for increasingly/decreasingly ordered selectionNT = 5;
NR = 4;
I = eye(NR,NR);
SNRs_dB = 0:2:20;
SNRs = 10.^(SNRs_dB/10);
N_SNR = length(SNRs);
Niter = 1000;
ns_antenna_select= [1 2 3 4];
N_ant = length(ns_antenna_select);
gss=['-ko';'-k^';'-kd';'-ks'];%% 主函数
for iant = 1:N_antfor isnr = 1:N_SNRSNR = SNRs(isnr)/ns_antenna_select(iant);sum_capacity = 0;for i = 1:Niterif method_sel==0sel_ant_indices=[];rem_ant_indices=1:NT;elsesel_ant_indices=1:NT;endH = (randn(NR,NT)+1j*randn(NR,NT))/sqrt(2);if method_sel == 0for n_antenna_select = 1:ns_antenna_select(iant)clear log_SH;for ii = 1:length(rem_ant_indices)H_sel = H(:,[sel_ant_indices rem_ant_indices(ii)]);log_SH(ii) = log2(real(det(I+SNR*H_sel*H_sel')));end[max_cap, max_index] = max(log_SH);sel_ant_index = rem_ant_indices(max_index);rem_ant_indices = [rem_ant_indices(1:max_index-1) rem_ant_indices(max_index+1:end)];sel_ant_indices = [sel_ant_indices sel_ant_index];endelsefor n_antenna_select = 1:NT-ns_antenna_select(iant)clear log_SH;for ii=1:length(sel_ant_indices)H_sel = H(:,[sel_ant_indices(1:ii-1) sel_ant_indices(ii+1:end)]);log_SH(ii) = log2(real(det(I+SNR*H_sel*H_sel')));end[max_cap, max_index] = max(log_SH);sel_ant_indices = [sel_ant_indices(1:max_index-1) sel_ant_indices(max_index+1:end)];endendsum_capacity = sum_capacity + max_cap;endsel_capacity(iant,isnr) = sum_capacity/Niter;endplot(SNRs_dB,sel_capacity(iant,:),gss(ns_antenna_select(iant),:), 'LineWidth',2); hold on;
end

%-----------------------OSTBC的天线选择---------------------%
%-----------------------author:lzx--------------------------%
%-----------------------date:2022年4月30日09点58分-----------------%
%% 参数设置
clear;
clf;
L_frame = 1000;
Niter = 400;
Nmod = 2;
M = 2^Nmod;N_T_Alamouti = 4;
NT = 2;
NR = 1;SNRs_dB = 0:2:20;
SNRs = 10.^(SNRs_dB/10);
N_SNR = length(SNRs);BERs = zeros(2, N_SNR);
%% 主函数
for isnr = 1:N_SNRSNR = SNRs(isnr);sigma = sqrt(NT/(2*SNR));n_biterror = 0;for iiter = 1:Niter% 生成数据frame_origin = randi([0,1],L_frame,NT*Nmod);% QPSK调制frame_mod=QPSKMod(frame_origin,L_frame, NT);frame_transmit = reshape(frame_mod, NT, 1, L_frame);frame_alamouti = [frame_transmit(1,1, :) -conj(frame_transmit(2,1, :));frame_transmit(2,1, :) conj(frame_transmit(1,1, :))];% 生成信道，SIMO有NR个信道Hiid = (randn(NR,N_T_Alamouti)+1j*randn(NR,N_T_Alamouti))/sqrt(2);noise = sigma*(randn(NR,NT)+1j*randn(NR,NT));for ih = 1:N_T_Alamoutifro_h(ih) = norm(Hiid(:,ih),'fro');end[val,index] = sort(fro_h,'descend');H_sel = Hiid(:,index([1 2]));for iframe = 1:L_framey(:,:,iframe) = H_sel*frame_alamouti(:,:,iframe)+noise;end% 接收H_sel_re = [conj(H_sel(1)) H_sel(2); conj(H_sel(2)) -H_sel(1)];y_re = [y(1,1,:); conj(y(1,2,:))];for iframe = 1:L_frameframe_re(:,:,iframe) = (H_sel_re* y_re(:,:,iframe))./(norm(H_sel_re,"fro")^2);endframe_pre_demod = reshape(frame_re, L_frame, NT);frame_demod = QPSKDemod(frame_pre_demod,L_frame,NT);% 计算误码率n_biterror_tmp = sum(sum(abs(frame_demod - frame_origin)))n_biterror = n_biterror + n_biterror_tmp;endBERs(isnr) = n_biterror/(Niter*L_frame*Nmod*NT);
end
%% 画图
BERs(2, :) = [0.177905000000000 0.131740000000000   0.0911950000000000  0.0580825000000000  0.0335800000000000  0.0178250000000000  0.00892500000000000 0.00408000000000000 0.00184000000000000 0.000750000000000000    0.000295000000000000];
semilogy(SNRs_dB,BERs(1, :),"-^", SNRs_dB,BERs(2, :),"-x");
axis([SNRs_dB([1 end]) 1e-6 1e0]);
xlabel('SNR[dB]'), ylabel('BER');
legend('Precoded Alamouti','No Precoded Alamouti');

3.一些思考

这里的代码就没有很认真的写了，主要来自于https://zhuanlan.zhihu.com/p/363970679，感谢陈老湿
天线选择是个挺简单的技术，首先你本身发射的要用到的天线个数就是小于总天线个数的

如何选择呢

最简单的就是做排列组合啊， $C_{n}^{M}$ ,每一个都算，算出来最大的就行了，但这样复杂度太高了，因此有了两个降低复杂度的算法

第一种是找出最好的，先找出最好的一根，然后吧剩下的和这一根组合，算（1，2）（1，3）（1，4）看看谁最大，然后一个个加

第二种是删除最差的，首先都选上，然后先删去一个，看看删去哪一个后剩下的最大就删去他，删完之后重复删，一直删到目标。

OSTBC选择技术

成对错差概率是完全看不懂啊

不过仿真非常简单

就是找出Fro范数最大的信道使用就行了，得到的结果也确实好一些

总结

这一章本以为能学到很多，其实并没有，我最想知道的CSI反馈流程没有体现出来

后面找机会再学习吧！