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

最新推荐文章于 2024-05-04 09:40:42 发布

举熙熙然回巧献计

最新推荐文章于 2024-05-04 09:40:42 发布

阅读量3.8k

点赞数 12

分类专栏： MIMO 文章标签： matlab 开发语言 dsp开发网络协议 ofdma

本文链接：https://blog.csdn.net/weixin_44863193/article/details/124549393

版权

MIMO 专栏收录该内容

7 篇文章 66 订阅

订阅专栏

初识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_SNR
    SNR = SNRs(isnr);
    sigma = sqrt(NT/(2*SNR));
    n_biterror = 0;
    for iiter = 1:N_iter
        frame_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_code
            arg_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_frame
            Rx(:,:,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_frame
            y(:,:,iframe) = HW_re* Rx_re(:,:,iframe)./norm(HW)^2;
        end
        frame_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;
    end
    BERs(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 ==8
    cloumn_index = [1];
    rotation_vector = [1,0];
elseif NT == 3 && M == 1 && L ==32
    cloumn_index=[1]; 
    rotation_vector=[1 26 28];
elseif NT == 4 && M == 2 && L ==32
    cloumn_index=[1 2]; 
    rotation_vector=[1 26 28];
elseif NT == 4 && M ==1 && L == 64
    cloumn_index=[1]; 
    rotation_vector=[1 8 61 45];
elseif NT == 4 && M ==2 && L == 64
    cloumn_index=[0 1]; 
    rotation_vector=[1 7 52 56];
elseif NT == 4 && M ==3 && L == 64
    cloumn_index=[0 2 3]; 
    rotation_vector=[1 8 61 45];
else
    error("没有对应的码本");
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-1
    code_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% ZF
        W_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% MMSE
        W_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 == 4
        W_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)))')));
    end
    for isnr = 1:N_SNR
        SNR = 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;
%                 end
                y(:,:,iframe) = (Hiid * W_pre * frame_transmit(:,:,iframe)+noise)/beta;
            end
            % 信号检测
            W = W_formula(Hiid, sigma, NT);
            for iframe = 1:L_frame
                x_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;
        end
        BERs(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:NT
    if n_ant>NT || n_ant<1
       error('n_ant must be between 1 and NT!');
     else
        Types = nchoosek([1:NT],n_ant);
        N_type = size(Types,1);
     end
   for isnr=1:N_SNR
      SNR = SNRs(isnr)/n_ant;  
%       rand('seed',1); randn('seed',1);  
      sum_capacity = 0;
      for iiter=1:N_iter
         H = (randn(NR,NT)+1j*randn(NR,NT))/sqrt(2);
         log_SH = zeros(1,N_type);
         for itype=1:N_type
            H_sel = H(:,Types(itype,:)); 
            log_SH(itype)=log2(real(det(I+SNR*H_sel*H_sel'))); % Eq.(12.22)
         end
         sum_capacity = sum_capacity + max(log_SH);
      end
      Capacitys(isnr) = sum_capacity/N_iter;
   end
   plot(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 selection

NT = 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_ant
    for isnr = 1:N_SNR
        SNR = SNRs(isnr)/ns_antenna_select(iant);
        sum_capacity = 0;
        for i = 1:Niter
            if method_sel==0
                sel_ant_indices=[];
                rem_ant_indices=1:NT;
            else
                sel_ant_indices=1:NT;
            end
            H = (randn(NR,NT)+1j*randn(NR,NT))/sqrt(2);
            if method_sel == 0
                for 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];
                end
            else
                for 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)];
                end
            end
            sum_capacity = sum_capacity + max_cap;
        end
        sel_capacity(iant,isnr) = sum_capacity/Niter;
    end
    plot(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_SNR
    SNR = 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_Alamouti
            fro_h(ih) = norm(Hiid(:,ih),'fro');
        end
        [val,index] = sort(fro_h,'descend');
        H_sel = Hiid(:,index([1 2]));
        for iframe = 1:L_frame
            y(:,:,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_frame
            frame_re(:,:,iframe) = (H_sel_re* y_re(:,:,iframe))./(norm(H_sel_re,"fro")^2);
        end
        frame_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;
    end
    BERs(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反馈流程没有体现出来

后面找机会再学习吧！

举熙熙然回巧献计

关注

12
点赞
踩
42

收藏

觉得还不错? 一键收藏
1
评论
初识MIMO(五)：CSI反馈及其仿真

初识MIMO(五)：CSI反馈及其仿真零代码地址https://github.com/liu-zongxi/MIMO_simulation一. 发射端的信道估计这章其实给了我隔靴搔痒的感觉，CSI的估计要尽可能快的完成，这里给出了这么几条路利用信道的互异性，这样可以在发射端估计量化CSI ，就是把CSI压缩后传输利用码本，码本表示CSI具体怎么做呢，书中没有给出，以后慢慢读论文吧！二. CSI 在接收端的利用技术——Alamouti预编码及其仿真1. 概念2.代码展示
复制链接

扫一扫