Returns the gain matrices of N stations with K receivers. function xsol = radio_interferometric_calibration(N, K) N >= K is always assumed. The example considers calibration of an array of N stations. We simulate a system with N stations, each having K receivers. For radio astronomy, K = 2. For a detailed exposition of the problem at hand, refer to the paper: "Radio interferometric calibration using a Riemannian manifold", Sarod Yatawatta, ICASSP, 2013. Available at http://dx.doi.org/10.1109/ICASSP.2013.6638382. The source of the signal is unpolarized (given by the matrix C). The measured data is the cross correlation of the signals at each receiver. So there will be N(N-1)/2 possible cross correlations. Noise with given SNR is added to the signal. The objective is to estimate the gains of each receiver (K x K) matrix, so the total size of the solutions is N x (K x K), which is written as an NK x K matrix. Note: each station gain matrix (KxK) can have a KxK unitary ambiguity, therefore we use the quotient manifold structure. The unitary ambiguity is common to all stations, so the solution obtained by optimization routine always has an unkown unitary matrix that makes the solution different from the true solution.
0001 function xsol = radio_interferometric_calibration(N, K) 0002 % Returns the gain matrices of N stations with K receivers. 0003 % 0004 % function xsol = radio_interferometric_calibration(N, K) 0005 % 0006 % N >= K is always assumed. 0007 % 0008 % The example considers calibration of an array of N stations. 0009 % We simulate a system with N stations, each having K receivers. 0010 % For radio astronomy, K = 2. 0011 % 0012 % For a detailed exposition of the problem at hand, refer to the paper: 0013 % "Radio interferometric calibration using a Riemannian manifold", 0014 % Sarod Yatawatta, ICASSP, 2013. 0015 % Available at http://dx.doi.org/10.1109/ICASSP.2013.6638382. 0016 % 0017 % The source of the signal is unpolarized (given by the matrix C). 0018 % The measured data is the cross correlation of the signals at each receiver. 0019 % So there will be N(N-1)/2 possible cross correlations. 0020 % Noise with given SNR is added to the signal. 0021 % 0022 % The objective is to estimate the gains of each receiver (K x K) matrix, 0023 % so the total size of the solutions is N x (K x K), which is written 0024 % as an NK x K matrix. 0025 % 0026 % Note: each station gain matrix (KxK) can have a KxK unitary ambiguity, 0027 % therefore we use the quotient manifold structure. The unitary ambiguity 0028 % is common to all stations, so the solution obtained by 0029 % optimization routine always has an unkown unitary matrix that makes the 0030 % solution different from the true solution. 0031 % 0032 0033 % This file is part of Manopt: www.manopt.org. 0034 % Original author: Sarod Yatawatta, June 29, 2015. 0035 % Contributors: Bamdev Mishra. 0036 % Change log: 0037 % 0038 % June 28, 2016 (BM): 0039 % Modified the egrad and ehess operations according to 0040 % the modified metric in the symfixedrankYYcomplexfactory file, 0041 % where a factor of 2 was removed from the metric. Accordingly, 0042 % a factor of 2 was added to egrad and ehess operations. 0043 % Aug 31, 2021 (XJ): 0044 % Added AD to compute the egrad and the ehess 0045 0046 % Generate some random data to test the function 0047 0048 if ~exist('N', 'var') || isempty(N) 0049 N = 10; 0050 end 0051 if ~exist('K', 'var') || isempty(K) 0052 K = 2; 0053 end 0054 0055 assert(N >= K, 'N must be larger than or equal to K.'); 0056 0057 % Baselines (pairs of correlations) 0058 B = N*(N-1)/2; 0059 0060 0061 0062 % Source coherence, at phase center 0063 C = eye(K); 0064 0065 % Random J (gains) of all stations 0066 J = 0.2*rand(K*N,K) + 1i*rand(K*N,K); 0067 0068 % Visibilities (cross correlations) 0069 V = zeros(K*B,K); 0070 0071 ck = 1; 0072 for ci = 1 : N -1 0073 for cj = ci + 1 : N 0074 % Compute cross correlation of each receiver pair. 0075 V(K*(ck-1)+1:K*ck,:) = J(K*(ci-1)+1:K*ci,:)*C*J(K*(cj-1)+1:K*cj,:)'; 0076 ck = ck + 1; 0077 end 0078 end 0079 0080 % Generate noise 0081 SNR = 10000;% inf; 0082 nn = randn(K*B,K)+1i*randn(K*B,K); 0083 noise_var = norm(V)^2/(norm(nn)^2*SNR); 0084 nn = nn*sqrt(noise_var); 0085 0086 % Add noise to signal 0087 V = V + nn; 0088 0089 0090 % Optimization part by creating the problem structure. 0091 % First, we use the manifold desctription. 0092 % Second, we define the problem cost, gradient and Hessian functions. 0093 0094 0095 % Manifold description 0096 % Note that the actual dimension is KN x K. 0097 problem.M = symfixedrankYYcomplexfactory(K*N, K); 0098 0099 0100 % Cost function 0101 problem.cost = @cost; 0102 function fval = cost(x) 0103 fval = 0.0; 0104 ck = 1; 0105 for p = 1 : N - 1 0106 for q = p + 1 : N 0107 res = V(K*(ck-1)+1:K*ck,:) - x(K*(p-1)+1:K*p,:)*C*x(K*(q-1)+1:K*q,:)'; % Residual 0108 fval = fval + real(res(:)'*res(:)); % Add norm of the residual. 0109 ck = ck + 1; 0110 end 0111 end 0112 end 0113 0114 % Euclidean gradient of the cost function. 0115 % Manopt automatically converts it to the Riemannian couterpart. 0116 % The code involves for-loops for readability, but could be vectorized 0117 % for improved speed. 0118 problem.egrad = @egrad; 0119 function grad = egrad(x) 0120 grad = zeros(K*N, K); 0121 ck = 1; 0122 for p = 1 : N - 1 0123 for q = p+1 : N 0124 res = 2*(V(K*(ck-1)+1:K*ck,:) - x(K*(p-1)+1:K*p,:)*C*x(K*(q-1)+1:K*q,:)'); % Residual 0125 grad(K*(p-1)+1:K*p,:) = grad(K*(p-1)+1:K*p,:) - res*x(K*(q-1)+1:K*q,:)*C'; 0126 grad(K*(q-1)+1:K*q,:) = grad(K*(q-1)+1:K*q,:) - res'*x(K*(p-1)+1:K*p,:)*C; 0127 ck = ck + 1; 0128 end 0129 end 0130 end 0131 0132 % Euclidean Hessian of the cost function along a search direction eta. 0133 % Manopt automatically converts it to the Riemannian couterpart. 0134 problem.ehess = @ehess; 0135 function hess = ehess(x, eta) 0136 hess = zeros(K*N, K); 0137 ck = 1; 0138 for p = 1 : N-1 0139 for q = p+1:N 0140 res = 2*(V(K*(ck-1)+1:K*ck,:) -x(K*(p-1)+1:K*p,:)*C*x(K*(q-1)+1:K*q,:)'); % Residual 0141 resdot = 2*(-x(K*(p-1)+1:K*p,:)*C*eta(K*(q-1)+1:K*q,:)' - eta(K*(p-1)+1:K*p,:)*C*x(K*(q-1)+1:K*q,:)'); % Residual derivative 0142 0143 hess(K*(p-1)+1:K*p,:) = hess(K*(p-1)+1:K*p,:) - (res*eta(K*(q-1)+1:K*q,:) + resdot*x(K*(q-1)+1:K*q,:))*C'; 0144 hess(K*(q-1)+1:K*q,:) = hess(K*(q-1)+1:K*q,:) - (res'*eta(K*(p-1)+1:K*p,:) + resdot'*x(K*(p-1)+1:K*p,:))*C; 0145 ck = ck + 1; 0146 end 0147 end 0148 end 0149 0150 0151 % An alternative way to compute the egrad and the ehess is to use 0152 % automatic differentiation provided in the deep learning toolbox. 0153 % Notice that the for loop in the cost function can make AD much 0154 % slower. Call manoptAD to prepare AD for the problem structure. 0155 % problem = manoptAD(problem); 0156 0157 0158 % Execute some checks on the derivatives for early debugging. 0159 % checkgradient(problem); 0160 % pause; 0161 % checkhessian(problem); 0162 % pause; 0163 0164 0165 % Solve. 0166 [xsol, xcost, info] = trustregions(problem); 0167 fprintf('Final cost: %g.\n', xcost); 0168 0169 0170 % Display some statistics. 0171 fs = 11; 0172 figure; 0173 semilogy([info.iter], [info.gradnorm], 'o-.','Color','blue', 'MarkerSize',6, 'LineWidth',1.1); 0174 ax1 = gca; 0175 set(ax1,'FontSize',fs); 0176 xlabel(ax1, 'Iteration #', 'FontSize',fs); 0177 ylabel(ax1, 'Gradient norm', 'FontSize',fs); 0178 title('Convergence of the trust-regions algorithm'); 0179 0180 % Make a plot of estimation error (only for K = 2). 0181 if K == 2 0182 % Find unitary ambiguity first by solving min ||J - xsol U||. 0183 % This has a closed-form solution. 0184 [u, ignore, v] = svd(xsol'*J); %#ok<ASGLU> 0185 0186 % Error in position 0187 E = J - xsol*u*v'; 0188 0189 % Normalize error 0190 E = E/norm(J); 0191 0192 % Plot 0193 figure; 0194 ax1 = subplot(1,2,1); 0195 quiver(real(J(:,1)), imag(J(:,1)),real(E(:,1)),imag(E(:,1))); 0196 hold all; 0197 scatter(real(J(:,1)), imag(J(:,1))); 0198 set(ax1,'FontSize',fs); 0199 xlabel('Real E_1'); 0200 ylabel('Imag E_1'); 0201 title('Position error 1st coordinate'); 0202 axis equal; 0203 ax2 = subplot(1,2,2); 0204 quiver(real(J(:,2)),imag(J(:,2)),real(E(:,2)),imag(E(:,2))); 0205 hold all; 0206 scatter(real(J(:,2)),imag(J(:,2))); 0207 set(ax2,'FontSize',fs); 0208 xlabel('Real E_2'); 0209 ylabel('Imag E_2'); 0210 title('Position error 2nd coordinate'); 0211 axis equal; 0212 end 0213 0214 end