mimo_capacity_water_filling.m

注水算法的系统的系统容量的仿真程序。对学习MIMO技术很有帮助。

% 程序-6.2 （2008.毕业设计-MIMO系统信道容量分析）
% 自适应功率分配的MIMO信道容量
% 发射端已知信道的状态信息时的信道容量 与 发射端未知信道的状态信息时的信道容量 的仿真对比
% 瑞利衰落信道矩阵
% MIMO_Capacity_water_filling.m

clear all;clc;

Nt = 3;											                           % 发射天线数目
Nr = 3;											                           % 接收天线数目

SampleNum = 10000;								                           % 蒙特卡洛仿真抽样数目
SNR_dB = 0:2:16;								                           % 信噪比（单位dB)

for nSNR = 1:length(SNR_dB)
	
	SNR_dB(nSNR)
	rho = 10^(SNR_dB(nSNR)/10);                                            % 信噪比的转换关系式
	
	for nSample = 1:SampleNum

		H = (randn(Nr,Nt)+sqrt(-1)*randn(Nr,Nt))/sqrt(2);                  % 瑞利衰落信道矩阵
		
		C_CSIR(nSNR,nSample) = log2(real(det(eye(Nt)+rho/Nt*H'*H)));
		
		EigCH = real(eig(H*H')).';                                         % 计算信道矩阵的特征值
		EigCH = EigCH(find(EigCH~=0));
		EQNoiseVar = 1./EigCH;
		
		[Pow, NoiseVarOut] = Water_Filling(rho, EQNoiseVar);               % 注水算法
		
		C_CSIT(nSNR,nSample) = sum(log2(real(1+Pow./NoiseVarOut)));
		
	end
	
	C_Unknown1(nSNR) = mean(C_CSIR(nSNR,:));
	C_Known1(nSNR)   = mean(C_CSIT(nSNR,:));
	
end

Nt = 5;											                           % 发射天线数目
Nr = 5;											                           % 接收天线数目

SampleNum = 10000;								                           % 蒙特卡洛仿真抽样数目
SNR_dB = 0:2:16;								                           % 信噪比（单位dB)

for nSNR = 1:length(SNR_dB)
	
	SNR_dB(nSNR)
	rho = 10^(SNR_dB(nSNR)/10);                                            % 信噪比的转换关系式
	
	for nSample = 1:SampleNum

		H = (randn(Nr,Nt)+sqrt(-1)*randn(Nr,Nt))/sqrt(2);                  % 瑞利衰落信道矩阵
		
		C_CSIR(nSNR,nSample) = log2(real(det(eye(Nt)+rho/Nt*H'*H)));
		
		EigCH = real(eig(H*H')).';                                         % 计算信道矩阵的特征值
		EigCH = EigCH(find(EigCH~=0));
		EQNoiseVar = 1./EigCH;
		
		[Pow, NoiseVarOut] = Water_Filling(rho, EQNoiseVar);               % 注水算法
		
		C_CSIT(nSNR,nSample) = sum(log2(real(1+Pow./NoiseVarOut)));
		
	end
	
	C_Unknown2(nSNR) = mean(C_CSIR(nSNR,:));
	C_Known2(nSNR)   = mean(C_CSIT(nSNR,:));
	
end

Nt = 7;											                           % 发射天线数目
Nr = 7;											                           % 接收天线数目

SampleNum = 10000;								                           % 蒙特卡洛仿真抽样数目
SNR_dB = 0:2:16;								                           % 信噪比（单位dB)

for nSNR = 1:length(SNR_dB)
	
	SNR_dB(nSNR)
	rho = 10^(SNR_dB(nSNR)/10);                                            % 信噪比的转换关系式
	
	for nSample = 1:SampleNum

		H = (randn(Nr,Nt)+sqrt(-1)*randn(Nr,Nt))/sqrt(2);                  % 瑞利衰落信道矩阵
		
		C_CSIR(nSNR,nSample) = log2(real(det(eye(Nt)+rho/Nt*H'*H)));
		
		EigCH = real(eig(H*H')).';                                         % 计算信道矩阵的特征值
		EigCH = EigCH(find(EigCH~=0));
		EQNoiseVar = 1./EigCH;
		
		[Pow, NoiseVarOut] = Water_Filling(rho, EQNoiseVar);               % 注水算法
		
		C_CSIT(nSNR,nSample) = sum(log2(real(1+Pow./NoiseVarOut)));
		
	end
	
	C_Unknown3(nSNR) = mean(C_CSIR(nSNR,:));
	C_Known3(nSNR)   = mean(C_CSIT(nSNR,:));
	
end

Nt = 9;											                           % 发射天线数目
Nr = 9;											                           % 接收天线数目

SampleNum = 10000;								                           % 蒙特卡洛仿真抽样数目
SNR_dB = 0:2:16;								                           % 信噪比（单位dB)

for nSNR = 1:length(SNR_dB)
	
	SNR_dB(nSNR)
	rho = 10^(SNR_dB(nSNR)/10);                                            % 信噪比的转换关系式
	
	for nSample = 1:SampleNum

		H = (randn(Nr,Nt)+sqrt(-1)*randn(Nr,Nt))/sqrt(2);                  % 瑞利衰落信道矩阵
		
		C_CSIR(nSNR,nSample) = log2(real(det(eye(Nt)+rho/Nt*H'*H)));
		
		EigCH = real(eig(H*H')).';                                         % 计算信道矩阵的特征值
		EigCH = EigCH(find(EigCH~=0));
		EQNoiseVar = 1./EigCH;
		
		[Pow, NoiseVarOut] = Water_Filling(rho, EQNoiseVar);               % 注水算法
		
		C_CSIT(nSNR,nSample) = sum(log2(real(1+Pow./NoiseVarOut)));
		
	end
	
	C_Unknown4(nSNR) = mean(C_CSIR(nSNR,:));
	C_Known4(nSNR)   = mean(C_CSIT(nSNR,:));
	
end

plot(SNR_dB, C_Unknown1, 'k-');                                             %开始仿真作图
hold on;
plot(SNR_dB, C_Known1, 'r-.');
hold on;
plot(SNR_dB, C_Unknown2, 'k-');                                             %开始仿真作图
hold on;
plot(SNR_dB, C_Known2, 'r-.');
hold on;
plot(SNR_dB, C_Unknown3, 'k-');                                             %开始仿真作图
hold on;
plot(SNR_dB, C_Known3, 'r-.');
hold on;
plot(SNR_dB, C_Unknown4, 'k-');                                             %开始仿真作图
hold on;
plot(SNR_dB, C_Known4, 'r-.');
hold on;

xlabel('信噪比 (dB)')
ylabel('遍历容量')
legend('未知信道状态信息','已知信道状态信息');
text(5,70,'从下往上依次为3X3,5X5,7X7,9X9天线的信道')