-
Notifications
You must be signed in to change notification settings - Fork 30
/
WD_solve_T.m
46 lines (40 loc) · 1.58 KB
/
WD_solve_T.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
function T_solve = WD_solve_T(X_T12,Y_T12,X_T13,Y_T13)
% 该函数用来求解失真矩阵 T
%
% 输入变量:
% 1)X_T12 是 S_T12 的特征向量矩阵;
% 2)Y_T12 是 M_T12 的特征向量矩阵;
% 3)X_T13 是 S_T13 的特征向量矩阵;
% 4)Y_T13 是 M_T13 的特征向量矩阵;
%
% 输出值:
% T_solve 是求解得到的发射失真矩阵 T ;
%
% 该程序的截止至:2016.05.16. 21:10
%%
% 计算 c2/c1 的数值
% 根据公式(161a),如下:
c2_c1 = ((X_T12(1,1)*X_T13(2,1) - X_T12(2,1)*X_T13(1,1))...
*(Y_T12(2,2)*Y_T13(1,1) - Y_T12(1,2)*Y_T13(2,1)))...
/((X_T12(2,2)*X_T13(1,1) - X_T12(1,2)*X_T13(2,1))...
*(Y_T12(1,1)*Y_T13(2,1) - Y_T12(2,1)*Y_T13(1,1)));
%{
% 根据公式(161b),如下:
c2_c1_2 = ((X_T12(1,1)*X_T13(2,2) - X_T12(2,1)*X_T13(1,2))...
*(Y_T12(2,2)*Y_T13(1,2) - Y_T12(1,2)*Y_T13(2,2)))...
/((X_T12(2,2)*X_T13(1,2) - X_T12(1,2)*X_T13(2,2))...
*(Y_T12(1,1)*Y_T13(2,2) - Y_T12(2,1)*Y_T13(1,2)));
disp('-------------------------------------------------------------------')
disp('两种不同计算公式计算 c2/c1 的数值,它们的差异如下,理论上应该一致');
disp('c2_c1 - c2_c1_2 =');
c2_c1 - c2_c1_2 % 两种不同计算公式的结果的差异,理论上应该一致
disp('考察是否一致(在计算精度以内),若一致,可继续后续操作');
% 理论上,(161a)和(161b)的结果是一致的,在上述仿真中也得到了验证;
% 因此,实际计算中只需任选一种进行后续计算即可。
%}
% 下面根据上述结果计算 c1 和 c2 的数值
c1 = det(Y_T12)*(X_T12(1,1)*Y_T12(2,2) - c2_c1*X_T12(1,2)*Y_T12(2,1))^(-1);
c2 = det(Y_T12)*((1/c2_c1)*X_T12(1,1)*Y_T12(2,2) - X_T12(1,2)*Y_T12(2,1))^(-1);
% 将 c1 和 c2 的值代入式(157),即可求得矩阵 T,如下:
T_solve = X_T12*[c1,0;0,c2]*(inv(Y_T12));
end