matlab代码--基于注水法的MIMO信道容量实现

今天接触一个简单的注水法程序，搞懂数学原理即可看懂代码。

1 注水法简介

详细原理可以参考：
MIMO的信道容量以及实现
大致理论就是利用拉格朗日乘子法，求解信道容量的最大化问题，得到的解形如往水池中注水的形式，最后根据公式敲代码即可。

2 代码及结果

信道可以替换多种：
如常用的复高斯信道，正交信道，对角信道，获得相应的特征值进行信道容量计算

1 正交信道
B=randn(Nt,Nt)+i*randn(Nt,Nt); %产生信道
h1=orth(B); %获得正交信道
lambda = sort(eig(h1*h1'),'descend')'; %信道降序排序

2 对角信道
h_v = 1/sqrt(2)*(randn(1,Nt)+1j*randn(1,Nt)); %产生信道矩阵中对角线的值
h1 = diag(h_v,0); %产生对角信道矩阵
lambda = sort(eig(h1*h1'),'descend')'; %信道降序排序
3 常用的iid信道
h1 = 1/sqrt(2)*(randn(Nr(1),Nt(1))+1j*randn(Nr(1),Nt(1)));%添加iid信道
lambda = sort(eig(h1*h1'),'descend')'; %信道降序排序

完整仿真代码如下：

clc;
clear all；
close all;
Nr = 6; %接收天线
Nt = 6; %发送天线
snrdb = -5:1:10; %信噪比变化趋势，单位瓦特
snr = 10.^(snrdb/10); % 单位dB
Niterr = 1000; % 蒙特卡罗仿真次数
c = zeros(length(snr),1);%变量用于存储
mu = zeros(1,length(snr));%变量用于存储
p_opt = zeros(length(snr),2);%变量用于存储
for q =1:Niterr
h1 = 1/sqrt(2)*(randn(Nr(1),Nt(1))+1j*randn(Nr(1),Nt(1)));%添加iid信道
lambda = sort(eig(h1*h1'),'descend')'; %信道降序排序
for j = 1:length(snr)
%% 使用注水法
r = rank(h1);
    niter=0;
    while 1
        mu(j) = 1/(r-niter)*(1+1/snr(j)*sum(1./lambda(1:(r-niter))));  %功率门限值
        p_opt(j,1:(r-niter)) = mu(j)- 1./(snr(j)*lambda(1:(r-niter))); %注水法公式
    if abs(p_opt(j,:))== p_opt(j,:)
        break
    else
        p_opt(j,end-niter:end) = 0;        
        niter = niter+1;
    end
    end
    c(j,1) = c(j,1)+ sum(log2(1+snr(j)*p_opt(j,:).*lambda(1:rank(h1))));%信道容量公式
end
end
grid on;
plot(snrdb,c(:,1)/Niterr,'r-*');
xlabel('SNR(dB)');
ylabel('Channel Capacity (bits/sec)');
hold on;


%%%%%%%
clc;
Nr = 12; %接收天线
Nt = 12; %发送天线
snrdb = -5:1:10; %信噪比变化趋势，单位瓦特
snr = 10.^(snrdb/10); % 单位dB
Niterr = 1000; % 蒙特卡罗仿真次数
c = zeros(length(snr),1);%变量用于存储
mu = zeros(1,length(snr));%变量用于存储
p_opt = zeros(length(snr),2);%变量用于存储
for q =1:Niterr
h1 = 1/sqrt(2)*(randn(Nr(1),Nt(1))+1j*randn(Nr(1),Nt(1)));%添加复信道
lambda = sort(eig(h1*h1'),'descend')'; %信道降序排序
for j = 1:length(snr)
%% 使用注水法
r = rank(h1);
    niter=0;
    while 1
        mu(j) = 1/(r-niter)*(1+1/snr(j)*sum(1./lambda(1:(r-niter))));
        p_opt(j,1:(r-niter)) = mu(j)- 1./(snr(j)*lambda(1:(r-niter))); %注水法公式
    if abs(p_opt(j,:))== p_opt(j,:)
        break
    else
        p_opt(j,end-niter:end) = 0;        
        niter = niter+1;
    end
    end
    c(j,1) = c(j,1)+ sum(log2(1+snr(j)*p_opt(j,:).*lambda(1:rank(h1))));
end
end
grid on;
plot(snrdb,c(:,1)/Niterr,'k-s');
xlabel('SNR(dB)');
ylabel('Channel Capacity (bits/sec)');
hold on;

%%
clc;
Nr = 24; %接收天线
Nt = 24; %发送天线
snrdb = -5:1:10; %信噪比变化趋势，单位瓦特
snr = 10.^(snrdb/10); % 单位dB
Niterr = 1000; % 蒙特卡罗仿真次数
c = zeros(length(snr),1);%变量用于存储
mu = zeros(1,length(snr));%变量用于存储
p_opt = zeros(length(snr),2);%变量用于存储
for q =1:Niterr
h1 = 1/sqrt(2)*(randn(Nr(1),Nt(1))+1j*randn(Nr(1),Nt(1)));%添加复信道
lambda = sort(eig(h1*h1'),'descend')'; %信道降序排序
for j = 1:length(snr)
%% 使用注水法
r = rank(h1);
    niter=0;
    while 1
        mu(j) = 1/(r-niter)*(1+1/snr(j)*sum(1./lambda(1:(r-niter))));
        p_opt(j,1:(r-niter)) = mu(j)- 1./(snr(j)*lambda(1:(r-niter))); %注水法公式
    if abs(p_opt(j,:))== p_opt(j,:)
        break
    else
        p_opt(j,end-niter:end) = 0;        
        niter = niter+1;
    end
    end
    c(j,1) = c(j,1)+ sum(log2(1+snr(j)*p_opt(j,:).*lambda(1:rank(h1))));
end
end
%%画图 
grid on;
plot(snrdb,c(:,1)/Niterr,'g-o');
xlabel('SNR(dB)');
ylabel('Channel Capacity (bits/sec)');
hold on;
xticks([-5:1:10])
legend("Nt=6 Nr=6 MIMO", "Nt=12 Nr=12 MIMO", "Nt=24 Nr=24 MIMO")