## PURPOSE

Returns the gain matrices of N stations with K receivers.

## DESCRIPTION

``` Returns the gain matrices of N stations with K receivers.

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.```

## CROSS-REFERENCE INFORMATION

This function calls:
• symfixedrankYYcomplexfactory Manifold of n x n complex Hermitian pos. semidefinite matrices of rank k.
• norm NORM Norm of a TT/MPS tensor.
• norm NORM Norm of a TT/MPS block-mu tensor.
• trustregions Riemannian trust-regions solver for optimization on manifolds.
This function is called by:

## SOURCE CODE

```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):
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.
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
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.
0156
0157
0158     % Execute some checks on the derivatives for early debugging.
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);
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```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005