My_Whitt_sub.m

function [T_solve,R_solve,A_solve] = My_Whitt_sub(S1,M1,S2,M2,S3,M3,Flag1,Flag2)
% 我的改进Whitt算法的 子函数
%
% 参考《多极化合成孔径雷达定标技术研究》席育孝（硕士学位论文）；
% 但关键的：
%   1）对判断条件进行了修改（我的一种改进,利用系统失真矩阵的串扰大小）;
%   2）但这只是Whitt改进算法的步骤一,不完全，步骤二需利用额外的角反射器数据.
%
% 输入变量：
%   1）S1 和 M1 是三面角反射器的散射矩阵和观测矩阵；
%   2）S2 和 M2 是0°二面角反射器的散射矩阵和观测矩阵；
%   3）S3 和 M3 是45°二面角反射器的散射矩阵和观测矩阵；
%   默认为以上顺序，最好不要改变；
%   4）Flag1 是标记1，表示人为将 Y_T13 的特征向量交换次序；
%   5）Flag2 是标记2，表示人为将 Y_R13 的特征向量交换次序；
%   注：标记1和标记2，数值为0表示不做人为操作，数值为1表示进行人为操作；
%
% 输出变量：
%   1）T_solve 是求解得到的发射失真矩阵T;
%   2）R_solve 是求解得到的接收失真矩阵R;
%   3）A_solve 是求解得到的绝对幅度因子A;
%
% 该程序截止至：2016.10.10. 16:42

%%
% ***********************************************************
%
%               使用 Whitt 算法进行极化定标
%
% ***********************************************************

%%
% ----------------------------------------
% 求解 发射失真矩阵 T
% ----------------------------------------
% disp('-------------------------------------------------------------------')
% disp('下面求解 发射失真矩阵 T');

M_T12 = (M1)\M2;
S_T12 = (S1)\S2;

M_T13 = (M1)\M3;
S_T13 = (S1)\S3;

% 对 M_T12 和 S_T12 进行相应的特征值和特征向量求解
% 求取：
%   1） S_T12 的特征值 LAMDA_S_T12 及特征向量矩阵 X_T12；
%   2） M_T12 的特征值 LAMDA_M_T12 及特征向量矩阵 Y_T12；
[X_T12,LAMDA_S_T12] = eig(S_T12);
[Y_T12,LAMDA_M_T12] = eig(M_T12);
% 注意：此时 X_T12 和 Y_T12 的次序并不一定是对应的；
%       我暂时假定它是对应的并进行后续计算；
%       真正的判断要等到失真矩阵求解出来后再进行；

% 同理，对 M_T13 和 S_T13 也进行相应的特征值和特征向量求解
% 求取：
%   1） S_T13 的特征值 LAMDA_S_T13 及特征向量矩阵 X_T13；
%   2） M_T13 的特征值 LAMDA_M_T13 及特征向量矩阵 Y_T13；
[X_T13,LAMDA_S_T13] = eig(S_T13);
[Y_T13,LAMDA_M_T13] = eig(M_T13);
% 注意：此时 X_T13 和 Y_T13 的次序并不一定是对应的；
%       我暂时假定它是对应的并进行后续计算；
%       真正的判断要等到失真矩阵求解出来后再进行；

% 假定上述特征向量是对应的，进行第一次矩阵 T 的求解：
T_solve = WD_solve_T(X_T12,Y_T12,X_T13,Y_T13);

% 接下来是关键的条件判断：
%   如果上面的假定（特征向量是对应的）是不成立的；
%   则需要经过调整后重新计算失真矩阵的数值；
% 判断条件 1
if abs(T_solve(1,2)) > abs(T_solve(1,1))...
        || abs(T_solve(1,2)) > abs(T_solve(2,2))...
        || abs(T_solve(2,1)) > abs(T_solve(1,1))...
        || abs(T_solve(2,1)) > abs(T_solve(2,2))
    % 交换 LAMDA_M_T12 的特征值次序
    new_LAMDA_M_T12 = zeros(2,2);
    new_LAMDA_M_T12(1,1) = LAMDA_M_T12(2,2);
    new_LAMDA_M_T12(2,2) = LAMDA_M_T12(1,1);
    clear LAMDA_M_T12;
    LAMDA_M_T12 = new_LAMDA_M_T12;
    clear new_LAMDA_M_T12;
    % 交换特征向量矩阵 Y_T12 的次序（第1列和第2列互换）
    new_Y_T12 = zeros(2,2);
    new_Y_T12 = [Y_T12(:,2),Y_T12(:,1)];
    clear Y_T12;
    Y_T12 = new_Y_T12;
    clear new_Y_T12;
    % 完成上述交换后重新计算失真矩阵 T 的数值
    T_solve = WD_solve_T(X_T12,Y_T12,X_T13,Y_T13);
end

% 在经过上述判断和相应的调整后，再做一次检查，如果还有问题，则报错 ！！！
if abs(T_solve(1,2)) > abs(T_solve(1,1))...
        || abs(T_solve(1,2)) > abs(T_solve(2,2))...
        || abs(T_solve(2,1)) > abs(T_solve(1,1))...
        || abs(T_solve(2,1)) > abs(T_solve(2,2))
    disp('！！！！！！！！！！！！！！！！！！！！！！！！！')
    disp('                  Error');
    disp('         发射失真矩阵 T 求解错误，请检查');
    disp('！！！！！！！！！！！！！！！！！！！！！！！！！')
    return;
end
% disp('-------------------------------------------------------------------')
% disp('发射失真矩阵 T 求解完成')
% disp('-------------------------------------------------------------------')

if Flag1 == 1
    % 数值为0表示不做人为操作，数值为1表示进行人为操作;
    % 人为将 Y_T13 的特征向量交换次序;
    Y_T13 = [Y_T13(:,2),Y_T13(:,1)];
    T_solve = WD_solve_T(X_T12,Y_T12,X_T13,Y_T13);
end


%%
% ----------------------------------------
% 求解 接收失真矩阵 R
% ----------------------------------------
% 求解方法与上面类似
% disp('-------------------------------------------------------------------')
% disp('下面求解 接收失真矩阵 R')

M_R12 = M2/(M1);
S_R12 = S2/(S1);

M_R13 = M3/(M1);
S_R13 = S3/(S1);

% 对 M_R12 和 S_R12 进行相应的特征值和特征向量求解
% 求取：
%   1） S_R12 的特征值 LAMDA_S_R12 及特征向量矩阵 X_R12；
%   2） M_R12 的特征值 LAMDA_M_R12 及特征向量矩阵 Y_R12；
[X_R12,LAMDA_S_R12] = eig(S_R12);
[Y_R12,LAMDA_M_R12] = eig(M_R12);
% 注意：此时 X_R12 和 Y_R12 的次序并不一定是对应的；
%       我暂时假定它是对应的并进行后续计算；
%       真正的判断要等到失真矩阵求解出来后再进行；

% 同理，对 M_R13 和 S_R13 也进行相应的特征值和特征向量求解
% 求取：
%   1） S_R13 的特征值 LAMDA_S_R13 及特征向量矩阵 X_R13；
%   2） M_R13 的特征值 LAMDA_M_R13 及特征向量矩阵 Y_R13；
[X_R13,LAMDA_S_R13] = eig(S_R13) ;
[Y_R13,LAMDA_M_R13] = eig(M_R13);
% 注意：此时 X_R13 和 Y_R13 的次序并不一定是对应的；
%       我暂时假定它是对应的并进行后续计算；
%       真正的判断要等到失真矩阵求解出来后再进行；

% 假定上述特征向量是对应的，进行第一次矩阵 R 的求解：
R_solve = WD_solve_R(X_R12,Y_R12,X_R13,Y_R13);

% 接下来是关键的条件判断：
%   如果上面的假定（特征向量是对应的）是不成立的；
%   则需要经过调整后重新计算失真矩阵的数值；
% 判断条件 1
if abs(R_solve(1,2)) > abs(R_solve(1,1))...
        || abs(R_solve(1,2)) > abs(R_solve(2,2))...
        || abs(R_solve(2,1)) > abs(R_solve(1,1))...
        || abs(R_solve(2,1)) > abs(R_solve(2,2))
    % 交换 LAMDA_M_T12 的特征值次序
    new_LAMDA_M_R12 = zeros(2,2);
    new_LAMDA_M_R12(1,1) = LAMDA_M_R12(2,2);
    new_LAMDA_M_R12(2,2) = LAMDA_M_R12(1,1);
    clear LAMDA_M_R12;
    LAMDA_M_R12 = new_LAMDA_M_R12;
    clear new_LAMDA_M_R12;
    % 交换特征向量矩阵 Y_R12 的次序（第1列和第2列互换）
    new_Y_R12 = zeros(2,2);
    new_Y_R12 = [Y_R12(:,2),Y_R12(:,1)];
    clear Y_R12;
    Y_R12 = new_Y_R12;
    clear new_Y_T12;
    % 完成上述交换后重新计算失真矩阵 R 的数值
    R_solve = WD_solve_R(X_R12,Y_R12,X_R13,Y_R13);
end

% 在经过上述判断和相应的调整后，再做一次检查，如果还有问题，则报错 ！！！
if abs(R_solve(1,2)) > abs(R_solve(1,1))...
        || abs(R_solve(1,2)) > abs(R_solve(2,2))...
        || abs(R_solve(2,1)) > abs(R_solve(1,1))...
        || abs(R_solve(2,1)) > abs(R_solve(2,2))
    disp('！！！！！！！！！！！！！！！！！！！！！！！！！')
    disp('                  Error');
    disp('         接收失真矩阵 R 求解错误，请检查');
    disp('！！！！！！！！！！！！！！！！！！！！！！！！！')
    return;
end
% disp('-------------------------------------------------------------------')
% disp('接收失真矩阵 R 求解完成')
% disp('-------------------------------------------------------------------')

if Flag2 == 1
    % 数值为0表示不做人为操作，数值为1表示进行人为操作;
    % 人为将 Y_R13 的特征向量交换次序;    
    Y_R13 = [Y_R13(:,2),Y_R13(:,1)];
    R_solve = WD_solve_R(X_R12,Y_R12,X_R13,Y_R13);
end


%%
% ----------------------------------------
% 求解绝对幅度因子 A
% ----------------------------------------
% 任选三个目标中的一个目标，由散射矩阵中的某元素来计算绝对幅度因子 A

% 选取标准：
%   选择理论散射矩阵最准确的定标目标，可以获得 A 的最好估计
%
% 因为三面角反射器的散射矩阵是和视线方向无关的，而二面角则对视线方向角敏感；
% 所以我认为选择三面角反射器会得到更准确的估计，下面就基于三面角反射器进行计算；

% 用计算得到的 R_solve 和 T_solve 以及三面角反射器的散射矩阵 S1 计算得到
% 除绝对相位和绝对幅度外的矩阵数值
Tmp = R_solve*S1*T_solve;
% 任选矩阵的第 mn 个元素进行计算；
% 这里我选择的是（1,1）
A_solve = abs(M1(1,1))/abs(Tmp(1,1));

% disp('-------------------------------------------------------------------')
% disp('绝对幅度因子 A 求解完成');
% disp('-------------------------------------------------------------------')
% disp('用“Whitt算法”进行极化定标参数求解到此结束');
% disp('-------------------------------------------------------------------')

end