WD_solve_R.m

function R_solve = WD_solve_R(X_R12,Y_R12,X_R13,Y_R13)
% 该函数用来求解失真矩阵 R 
%
% 输入变量：
%   1）X_R12 是 S_R12 的特征向量矩阵；
%   2）Y_R12 是 M_R12 的特征向量矩阵；
%   3）X_R13 是 S_R13 的特征向量矩阵；
%   4）Y_R13 是 M_R13 的特征向量矩阵；
%
% 输出值：
%   R_solve 是求解得到的接收失真矩阵 R ;
%
% 该程序的截止至：2016.05.16. 21:14

%%
% 计算 d2/d1 的数值
% 根据公式（171a），如下：
d2_d1 = ((Y_R12(1,1)*Y_R13(2,1) - Y_R12(2,1)*Y_R13(1,1))...
    *(X_R12(2,2)*X_R13(1,1) - X_R12(1,2)*X_R13(2,1)))...
    /((Y_R12(2,2)*Y_R13(1,1) - Y_R12(1,2)*Y_R13(2,1))...
    *(X_R12(1,1)*X_R13(2,1) - X_R12(2,1)*X_R13(1,1)));

%{
% 根据公式（171b），如下：
d2_d1_2 = ((Y_R12(1,1)*Y_R13(2,2) - Y_R12(2,1)*Y_R13(1,2))...
    *(X_R12(2,2)*X_R13(1,2) - X_R12(1,2)*X_R13(2,2)))...
    /((Y_R12(2,2)*Y_R13(1,2) - Y_R12(1,2)*Y_R13(2,2))...
    *(X_R12(1,1)*X_R13(2,2) - X_R12(2,1)*X_R13(1,2)));

disp('-------------------------------------------------------------------')
disp('两种不同计算公式计算 d2/d1 的数值，它们的差异如下，理论上应该一致');
disp('d2_d1 - d2_d1_2 =');
d2_d1 - d2_d1_2     % 两种不同计算公式的结果的差异，理论上应该一致
disp('考察是否一致（在计算精度以内），若一致，可继续后续操作');
% 理论上，（171a）和（171b）的结果是一致的，在上述仿真中也得到了验证；
% 因此，实际计算中只需任选一种进行后续计算即可。
%}

% 下面根据上述结果计算 d1 和 d2 的数值
d1 = det(X_R12)*(Y_R12(1,1)*X_R12(2,2) - d2_d1*Y_R12(1,2)*X_R12(2,1))^(-1);
d2 = det(X_R12)*((1/d2_d1)*Y_R12(1,1)*X_R12(2,2) - Y_R12(1,2)*X_R12(2,1))^(-1);

% 将 d1 和 d2 的值代入式（168），即可求得矩阵 R，如下：
R_solve = Y_R12*[d1,0;0,d2]*(inv(X_R12));

end