Home > examples > radio_interferometric_calibration.m

radio_interferometric_calibration

PURPOSE ^

Returns the gain matrices of N stations with K receivers.

SYNOPSIS ^

function xsol = radio_interferometric_calibration(N, K)

DESCRIPTION ^

 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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

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     
0044     % Generate some random data to test the function
0045     
0046     if ~exist('N', 'var') || isempty(N)
0047         N = 10; 
0048     end
0049     if ~exist('K', 'var') || isempty(K)
0050         K = 2; 
0051     end
0052     
0053     assert(N >= K, 'N must be larger than or equal to K.');
0054     
0055     % Baselines (pairs of correlations)
0056     B = N*(N-1)/2;
0057     
0058     
0059     
0060     % Source coherence, at phase center
0061     C = eye(K);
0062     
0063     % Random J (gains) of all stations
0064     J = 0.2*rand(K*N,K) + 1i*rand(K*N,K);
0065  
0066     % Visibilities (cross correlations)
0067     V = zeros(K*B,K);
0068     
0069     ck = 1;
0070     for ci = 1 : N -1,
0071         for cj = ci + 1 : N,
0072             % Compute cross correlation of each receiver pair.
0073             V(K*(ck-1)+1:K*ck,:) = J(K*(ci-1)+1:K*ci,:)*C*J(K*(cj-1)+1:K*cj,:)';
0074             ck = ck + 1;
0075         end
0076     end
0077     
0078     % Generate noise
0079     SNR = 10000;% inf;
0080     nn = randn(K*B,K)+1i*randn(K*B,K);
0081     noise_var = norm(V)^2/(norm(nn)^2*SNR);
0082     nn = nn*sqrt(noise_var);
0083     
0084     % Add noise to signal
0085     V = V + nn;
0086     
0087     
0088     % Optimization part by creating the problem structure.
0089     % First, we use the manifold desctription.
0090     % Second, we define the problem cost, gradient and Hessian functions.
0091    
0092     
0093     % Manifold description
0094     % Note that the actual dimension is KN x K.
0095     problem.M = symfixedrankYYcomplexfactory(K*N, K);
0096     
0097     
0098     % Cost function
0099     problem.cost = @cost;
0100     function fval = cost(x)
0101         fval = 0.0;
0102         ck = 1;
0103         for p = 1 : N - 1,
0104             for q = p + 1 : N,
0105                 res = V(K*(ck-1)+1:K*ck,:) - x(K*(p-1)+1:K*p,:)*C*x(K*(q-1)+1:K*q,:)'; % Residual
0106                 fval = fval + real(res(:)'*res(:)); % Add norm of the residual.
0107                 ck = ck + 1;
0108             end
0109         end
0110     end
0111     
0112     % Euclidean gradient of the cost function.
0113     % Manopt automatically converts it to the Riemannian couterpart.
0114     % The code involves for-loops for readability, but could be vectorized
0115     % for improved speed.
0116     problem.egrad = @egrad;
0117     function grad = egrad(x)
0118         grad = zeros(K*N, K);
0119         ck = 1;
0120         for p = 1 : N - 1,
0121             for q = p+1 : N,
0122                 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
0123                 grad(K*(p-1)+1:K*p,:) = grad(K*(p-1)+1:K*p,:) - res*x(K*(q-1)+1:K*q,:)*C';
0124                 grad(K*(q-1)+1:K*q,:) = grad(K*(q-1)+1:K*q,:) - res'*x(K*(p-1)+1:K*p,:)*C;
0125                 ck = ck + 1;
0126             end
0127         end
0128     end
0129     
0130     % Euclidean Hessian of the cost function along a search direction eta.
0131     % Manopt automatically converts it to the Riemannian couterpart.
0132     problem.ehess = @ehess;
0133     function hess = ehess(x, eta)
0134         hess = zeros(K*N, K);
0135         ck = 1;
0136         for p = 1 : N-1,
0137             for q = p+1:N,
0138                 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
0139                 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
0140                 
0141                 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';
0142                 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;
0143                 ck = ck + 1;
0144             end
0145         end
0146     end
0147     
0148     
0149     
0150     % Execute some checks on the derivatives for early debugging.
0151     % checkgradient(problem);
0152     % pause;
0153     % checkhessian(problem);
0154     % pause;
0155     
0156     
0157     % Solve.
0158     [xsol,  xcost,  info] = trustregions(problem); 
0159     fprintf('Final cost: %g.\n', xcost);
0160     
0161     
0162     % Display some statistics.
0163     fs = 11;
0164     figure;
0165     semilogy([info.iter], [info.gradnorm], 'o-.','Color','blue', 'MarkerSize',6, 'LineWidth',1.1);
0166     ax1 = gca;
0167     set(ax1,'FontSize',fs);
0168     xlabel(ax1, 'Iteration #', 'FontSize',fs);
0169     ylabel(ax1, 'Gradient norm', 'FontSize',fs);
0170     title('Convergence of the trust-regions algorithm');
0171 
0172     % Make a plot of estimation error (only for K = 2).
0173     if K == 2,
0174         % Find unitary ambiguity first by solving min ||J - xsol U||.
0175         % This has a closed-form solution.
0176         [u, ignore, v] = svd(xsol'*J); %#ok<ASGLU>
0177 
0178         % Error in position
0179         E = J - xsol*u*v'; 
0180 
0181         % Normalize error
0182         E = E/norm(J);
0183 
0184         % Plot
0185         figure;
0186         ax1 = subplot(1,2,1);
0187         quiver(real(J(:,1)), imag(J(:,1)),real(E(:,1)),imag(E(:,1)));
0188         hold all;
0189         scatter(real(J(:,1)), imag(J(:,1)));
0190         set(ax1,'FontSize',fs);
0191         xlabel('Real E_1');
0192         ylabel('Imag E_1');
0193         title('Position error 1st coordinate'); 
0194         axis equal;
0195         ax2 = subplot(1,2,2);
0196         quiver(real(J(:,2)),imag(J(:,2)),real(E(:,2)),imag(E(:,2)));
0197         hold all;
0198         scatter(real(J(:,2)),imag(J(:,2)));
0199         set(ax2,'FontSize',fs);
0200         xlabel('Real E_2');
0201         ylabel('Imag E_2');
0202         title('Position error 2nd coordinate'); 
0203         axis equal;
0204     end
0205     
0206 end

Generated on Sat 12-Nov-2016 14:11:22 by m2html © 2005