目錄

一. 介紹

二. MATLAB代碼

三. 運行結果與分析

---

一. 介紹

本篇將在MATLAB的仿真環(huán)境中對比MIMO幾種常見的信道估計方法的性能。

有關MIMO的介紹可看轉至此篇博客：

MIMO系統(tǒng)模型構建_嘮嗑！的博客-CSDN博客

在所有無線通信中，信號通過信道會出現(xiàn)失真，或者會添加各種噪聲。正確解碼接收到的信號就需要消除信道施加的失真和噪聲。為了弄清信道的特性，就需要信道估計。

信道估計有很多不同的方法，但是通用的流程可概括如下：

1. 設置一個數學模型，利用信道矩陣搭建起發(fā)射信號和接收信號之間的關系；
2. 發(fā)射已知信號（通常稱為參考信號或導頻信號）并檢測接收到的信號；
3. 通過對比發(fā)送信號和接收信號，確定信道矩陣中的每個元素。

二. MATLAB代碼

一共有四個代碼，包含三個函數代碼和一個主運行代碼。

主運行代碼用來產生最后的圖像

（1）main.m文件代碼

%由于信道數據隨機產生，每次運行出的圖像可能有略微差異%初始化
close all;
clear all;%%設定仿真參數rng('shuffle'); %產生隨機化種子，也可以使用另一函數randn('state',sum(100*clock));%設定蒙特卡洛仿真的數目
nbrOfMonteCarloRealizations = 1000;nbrOfCouplingMatrices = 50; %相關矩陣數目Nt = 8; %發(fā)射天線的數量，訓練序列的長度
Nr = 4; %接收天線的數量%訓練的總功率
totalTrainingPower_dB = 0:1:20; %單位為dB
totalTrainingPower = 10.^(totalTrainingPower_dB/10); %轉為線性范圍%最優(yōu)化算法
option = optimset('Display','off','TolFun',1e-7,'TolCon',1e-7,'Algorithm','interior-point');%比較不同的信道估計算法 
%實用蒙特卡洛仿真法
average_MSE_MMSE_estimator_optimal = zeros(length(totalTrainingPower),nbrOfCouplingMatrices,2); %最優(yōu)訓練下的MMSE估計法
average_MSE_MMSE_estimator_heuristic = zeros(length(totalTrainingPower),nbrOfCouplingMatrices,2); %啟發(fā)訓練下的MMSE估計法
average_MSE_MVU_estimator = zeros(length(totalTrainingPower),nbrOfCouplingMatrices,2); %最優(yōu)訓練下的MVU估計法
average_MSE_onesided_estimator = zeros(length(totalTrainingPower),nbrOfCouplingMatrices,2); %單邊線性估計法 
average_MSE_twosided_estimator = zeros(length(totalTrainingPower),nbrOfCouplingMatrices,2); %雙邊線性估計法%隨機信道統(tǒng)計量下的迭代
for statisticsIndex = 1:nbrOfCouplingMatrices%產生Weichselberger模型下的耦合矩陣V%元素均來自卡方分布（自由度為2）V = abs(randn(Nr,Nt)+1i*randn(Nr,Nt)).^2;V = Nt*Nr*V/sum(V(:)); %將矩陣Frobenius范數設為 Nt x Nr.%計算耦合矩陣的協(xié)方差矩陣R = diag(V(:));R_T = diag(sum(V,1)); %在Weichselberger模型下，計算發(fā)射端的協(xié)方差矩陣R_R = diag(sum(V,2)); %在Weichselberger模型下，計算接收端的協(xié)方差矩陣%使用MATLAB內置自帶的優(yōu)化算法，計算MMSE估計法下最優(yōu)的訓練功率分配trainingpower_MMSE_optimal = zeros(Nt,length(totalTrainingPower)); %每個訓練序列的功率分配向量for k = 1:length(totalTrainingPower) %遍歷每個訓練序列的功率分配trainingpower_initial = totalTrainingPower(k)*ones(Nt,1)/Nt; %初始設定功率均相等%使用fmincon函數來最優(yōu)化功率分配%最小化MSE，所有功率均非負trainingpower_MMSE_optimal(:,k) = fmincon(@(q) functionMSEmatrix(R,q,Nr),trainingpower_initial,ones(1,Nt),totalTrainingPower(k),[],[],zeros(Nt,1),totalTrainingPower(k)*ones(Nt,1),[],option);end%計算功率分配[eigenvalues_sorted,permutationorder] = sort(diag(R_T),'descend'); %計算和整理特征值[~,inversePermutation] = sort(permutationorder); %記錄特征值的orderq_MMSE_heuristic = zeros(Nt,length(totalTrainingPower));for k = 1:length(totalTrainingPower) %遍歷每個訓練功率alpha_candidates = (totalTrainingPower(k)+cumsum(1./eigenvalues_sorted(1:Nt,1)))./(1:Nt)'; %計算拉格朗日乘子的不同值optimalIndex = find(alpha_candidates-1./eigenvalues_sorted(1:Nt,1)>0 & alpha_candidates-[1./eigenvalues_sorted(2:end,1); Inf]<0); %找到拉格朗日乘子的αq_MMSE_heuristic(:,k) = max([alpha_candidates(optimalIndex)-1./eigenvalues_sorted(1:Nt,1) zeros(Nt,1)],[],2); %使用最優(yōu)的α計算功率分配endq_MMSE_heuristic = q_MMSE_heuristic(inversePermutation,:); %通過重新整理特征值來確定最終的功率分配%計算均勻功率分配q_uniform = (ones(Nt,1)/Nt)*totalTrainingPower;%蒙特卡洛仿真初始化vecH_realizations = sqrtm(R)*( randn(Nt*Nr,nbrOfMonteCarloRealizations)+1i*randn(Nt*Nr,nbrOfMonteCarloRealizations) ) / sqrt(2); %以向量的形式產生信道 vecN_realizations = ( randn(Nt*Nr,nbrOfMonteCarloRealizations)+1i*randn(Nt*Nr,nbrOfMonteCarloRealizations) ) / sqrt(2); %以向量的形式產生噪聲%對于每種估計方法計算MSE和訓練功率for k = 1:length(totalTrainingPower)%MMSE估計法：最優(yōu)訓練功率分配P_tilde = kron(diag(sqrt(trainingpower_MMSE_optimal(:,k))),eye(Nr)); %計算有效功率矩陣average_MSE_MMSE_estimator_optimal(k,statisticsIndex,1) = trace(R - (R*P_tilde'/(P_tilde*R*P_tilde' + eye(length(R))))*P_tilde*R); %計算MSEH_hat = (R*P_tilde'/(P_tilde*R*P_tilde'+eye(length(R)))) * (P_tilde*vecH_realizations+vecN_realizations); %使用蒙特卡洛仿真來計算該估計average_MSE_MMSE_estimator_optimal(k,statisticsIndex,2) = mean( sum(abs(vecH_realizations - H_hat).^2,1) ); %使用蒙特卡洛仿真來計算MSE%MMSE估計法：啟發(fā)式訓練功率分配MMSE P_tilde = kron(diag(sqrt(q_MMSE_heuristic(:,k))),eye(Nr));  %計算有效訓練矩陣average_MSE_MMSE_estimator_heuristic(k,statisticsIndex,1) = trace(R - (R*P_tilde'/(P_tilde*R*P_tilde' + eye(length(R))))*P_tilde*R); %計算MSEH_hat = (R*P_tilde'/(P_tilde*R*P_tilde'+eye(length(R)))) * (P_tilde*vecH_realizations + vecN_realizations); %使用蒙特卡洛仿真來計算該估計average_MSE_MMSE_estimator_heuristic(k,statisticsIndex,2) = mean( sum(abs(vecH_realizations - H_hat).^2,1) ); %使用蒙特卡洛仿真來計算MSE%MVY估計法: 最優(yōu)均勻訓練功率分配P_training = diag(sqrt(q_uniform(:,k))); %均勻功率分配P_tilde = kron(transpose(P_training),eye(Nr));  %計算有效訓練矩陣P_tilde_pseudoInverse = kron((P_training'/(P_training*P_training'))',eye(Nr)); %計算有效訓練矩陣的偽逆average_MSE_MVU_estimator(k,statisticsIndex,1) = Nt^2*Nr/totalTrainingPower(k); %計算MSEH_hat = P_tilde_pseudoInverse'*(P_tilde*vecH_realizations + vecN_realizations); %使用蒙特卡洛仿真來計算該估計average_MSE_MVU_estimator(k,statisticsIndex,2) = mean( sum(abs(vecH_realizations - H_hat).^2,1) ); %使用蒙特卡洛仿真來計算MSE%One-sided linear 估計法: 最優(yōu)訓練功率分配又被稱為 "LMMSE 估計法" P_training = diag(sqrt(q_MMSE_heuristic(:,k))); %使用最優(yōu)功率分配來計算訓練矩陣 P_tilde = kron(P_training,eye(Nr)); %計算有效訓練矩陣average_MSE_onesided_estimator(k,statisticsIndex,1) = trace(inv(inv(R_T)+P_training*P_training'/Nr)); %計算MSEAo = (P_training'*R_T*P_training + Nr*eye(Nt))\P_training'*R_T; %計算one-sided linear估計法中的矩陣A0 H_hat = kron(transpose(Ao),eye(Nr))*(P_tilde*vecH_realizations + vecN_realizations); %使用蒙特卡洛仿真來計算該估計average_MSE_onesided_estimator(k,statisticsIndex,2) = mean( sum(abs(vecH_realizations - H_hat).^2,1) );  %使用蒙特卡洛仿真來計算MS%Two-sided linear 估計法: 最優(yōu)訓練功率分配P_training = diag(sqrt(q_uniform(:,k))); %計算訓練矩陣，均勻功率分配P_tilde = kron(P_training,eye(Nr)); %計算有效訓練矩陣R_calE = sum(1./q_uniform(:,k))*eye(Nr); %計算協(xié)方差矩陣average_MSE_twosided_estimator(k,statisticsIndex,1) = trace(R_R-(R_R/(R_R+R_calE))*R_R); %計算MSEC1 = inv(P_training); %計算矩陣C1C2bar = R_R/(R_R+R_calE); %計算C2bar矩陣H_hat = kron(transpose(C1),C2bar)*(P_tilde*vecH_realizations + vecN_realizations);average_MSE_twosided_estimator(k,statisticsIndex,2) = mean( sum(abs(vecH_realizations - H_hat).^2,1) ); %使用蒙特卡洛仿真來計算MSendend%挑選訓練功率的子集
subset = linspace(1,length(totalTrainingPower_dB),5);normalizationFactor = Nt*Nr; %設定MSE標準化因子為trace(R), 標準化MSE為從0到1.%使用理論MSE公式畫圖
figure(1); hold on; box on;plot(totalTrainingPower_dB,mean(average_MSE_MVU_estimator(:,:,1),2)/normalizationFactor,'b:','LineWidth',2);plot(totalTrainingPower_dB,mean(average_MSE_twosided_estimator(:,:,1),2)/normalizationFactor,'k-.','LineWidth',1);
plot(totalTrainingPower_dB,mean(average_MSE_onesided_estimator(:,:,1),2)/normalizationFactor,'r-','LineWidth',1);plot(totalTrainingPower_dB(subset(1)),mean(average_MSE_MMSE_estimator_heuristic(subset(1),:,1),2)/normalizationFactor,'b+-.','LineWidth',1);
plot(totalTrainingPower_dB(subset(1)),mean(average_MSE_MMSE_estimator_optimal(subset(1),:,1),2)/normalizationFactor,'ko-','LineWidth',1);legend('MVU, optimal','Two-sided linear, optimal','One-sided linear, optimal','MMSE, heuristic','MMSE, optimal','Location','SouthWest')plot(totalTrainingPower_dB,mean(average_MSE_MMSE_estimator_heuristic(:,:,1),2)/normalizationFactor,'b-.','LineWidth',1);
plot(totalTrainingPower_dB,mean(average_MSE_MMSE_estimator_optimal(:,:,1),2)/normalizationFactor,'k-','LineWidth',1);
plot(totalTrainingPower_dB(subset),mean(average_MSE_MMSE_estimator_heuristic(subset,:,1),2)/normalizationFactor,'b+','LineWidth',1);
plot(totalTrainingPower_dB(subset),mean(average_MSE_MMSE_estimator_optimal(subset,:,1),2)/normalizationFactor,'ko','LineWidth',1);set(gca,'YScale','Log'); %縱軸為log范圍
xlabel('Total Training Power (dB)');
ylabel('Average Normalized MSE');
axis([0 totalTrainingPower_dB(end) 0.05 1]);title('Results based on theoretical formulas');%使用蒙特卡洛仿真畫理論運算圖
figure(2); hold on; box on;plot(totalTrainingPower_dB,mean(average_MSE_MVU_estimator(:,:,2),2)/normalizationFactor,'b:','LineWidth',2);
plot(totalTrainingPower_dB,mean(average_MSE_twosided_estimator(:,:,2),2)/normalizationFactor,'k-.','LineWidth',1);
plot(totalTrainingPower_dB,mean(average_MSE_onesided_estimator(:,:,2),2)/normalizationFactor,'r-','LineWidth',1);
plot(totalTrainingPower_dB(subset(1)),mean(average_MSE_MMSE_estimator_heuristic(subset(1),:,2),2)/normalizationFactor,'b+-.','LineWidth',1);
plot(totalTrainingPower_dB(subset(1)),mean(average_MSE_MMSE_estimator_optimal(subset(1),:,2),2)/normalizationFactor,'ko-','LineWidth',1);legend('MVU, optimal','Two-sided linear, optimal','One-sided linear, optimal','MMSE, heuristic','MMSE, optimal','Location','SouthWest')plot(totalTrainingPower_dB,mean(average_MSE_MMSE_estimator_heuristic(:,:,2),2)/normalizationFactor,'b-.','LineWidth',1);
plot(totalTrainingPower_dB,mean(average_MSE_MMSE_estimator_optimal(:,:,2),2)/normalizationFactor,'k-','LineWidth',1);
plot(totalTrainingPower_dB(subset),mean(average_MSE_MMSE_estimator_heuristic(subset,:,2),2)/normalizationFactor,'b+','LineWidth',1);
plot(totalTrainingPower_dB(subset),mean(average_MSE_MMSE_estimator_optimal(subset,:,2),2)/normalizationFactor,'ko','LineWidth',1);set(gca,'YScale','Log'); %縱軸為log范圍
xlabel('Total Training Power (dB)');
ylabel('Average Normalized MSE');
axis([0 totalTrainingPower_dB(end) 0.05 1]);title('Results based on Monte-Carlo simulations');

包含每行具體代碼的解釋

（2）三個函數文件

function [deviation,powerAllocation]=functionLagrangeMultiplier(eigenvaluesTransmitter,totalPower,k,alpha)
%Compute the MSE for estimation of the squared Frobenius norm of the
%channel matrix for a given training power allocation. 
%INPUT:
%eigenvaluesTransmitter = Vector with the active eigenvalues at the
%                         transmitter side
%totalPower             = Total power of the training sequence
%k                      = Vector with k parameter values 
%alpha                  = Langrange multiplier value
%
%OUTPUT:
%deviation              = Difference between available power and used power
%powerAllocation        = Training power allocation 
%Compute power allocation 
powerAllocation = sqrt(8*(1./alpha(:))*eigenvaluesTransmitter'/3).*cos(repmat((-1).^k*pi/3,[length(alpha) 1])-atan(sqrt(8*(1./alpha(:))*(eigenvaluesTransmitter.^3)'/27-1))/3)-repmat(1./eigenvaluesTransmitter',[length(alpha) 1]);%Deviation between total available power and the power that is used
deviation = abs(totalPower-sum(powerAllocation,2));

function MSE = functionMSEmatrix(R_diag,q_powerallocation,B)
%Compute the MSE for estimation of the channel matrix for a given training
%power allocation. 
%INPUT:
%R_diag            = Nt Nr x Nt Nr diagonal covariance matrix
%q_powerallocation = Nt x 1 vector with training power allocation
%B                 = Length of the training sequence.
%
%OUTPUT:
%MSE               = Mean Squared Error for estimation of the channel matrixP_tilde = kron(diag(sqrt(q_powerallocation)),eye(B));MSE = trace(R_diag - R_diag*(P_tilde'/(P_tilde*R_diag*P_tilde'+eye(length(R_diag))))*P_tilde*R_diag);

function MSE = functionMSEnorm(eigenvaluesTransmitter,eigenvaluesReceiver,powerAllocation)
%Compute the MSE for estimation of the squared Frobenius norm of the
%channel matrix for a given training power allocation. 
%INPUT:
%eigenvaluesTransmitter = Nt x 1 vector with eigenvalues at the
%                         transmitter side
%eigenvaluesReceiver    = Nr x 1 vector with eigenvalues at the
%                         receiver side
%powerAllocation        = Nt x 1 vector with training power allocation
%
%OUTPUT:
%MSE               = Mean Squared Error for estimation of the squared normMSE = sum(sum(((eigenvaluesTransmitter*eigenvaluesReceiver').^2 + 2*(powerAllocation.*eigenvaluesTransmitter.^3)*(eigenvaluesReceiver').^3)./(1+(powerAllocation.*eigenvaluesTransmitter)*eigenvaluesReceiver').^2));

注意：

• 此MIMO發(fā)射天線為8，接收天線為4；
• 三個函數文件的命名需要與函數保持一致；
• 先運行函數文件，再運行main主文件；
• 函數文件出現(xiàn)變量數報錯是正?，F(xiàn)象；
• 運行出來有兩個圖，選擇任意一個圖即可。

三. 運行結果與分析

分析：

1. 橫向看，當訓練功率增加時，均方誤差（MSE）在減小，符合信道估計的基本邏輯；
2. 縱向對比，MMSE（optimal）》MMSE(heuristic)》one-sided linear》two-side linear>MVU

》代表左邊優(yōu)于右邊，每一個位置代表一種信道估計方法