Added experiments for the sensor location error paper.

.gitignore

@@ -1,3 +1,4 @@
 *.asv
 *.fig
+.vscode
 experimental/

estimator/mle_sto_uc_1d.m

@@ -84,6 +84,7 @@
     x0(n + 1:end) = z;
 else
     x0 = x0(:);
+    % Unify to radians.
     switch (unit)
         case 'degree'
             x0(1:n) = deg2rad(x0(1:n));

examples/experiments/coarrays_music_crb/README.md

@@ -0,0 +1,27 @@
+## About
+
+This folder contains scripts for analyzing the performance of SS-MUSIC and DA-MUSIC for various sparse linear arrays. You can use these scripts to produce figures similar to those in my following papers:
+
+```
+@article{wang_coarrays_2017,
+    title = {Coarrays, {MUSIC}, and the {Cram\'{e}r}-{Rao} Bound},
+    volume = {65},
+    issn = {1053-587X},
+    doi = {10.1109/TSP.2016.2626255},
+    number = {4},
+    journal = {IEEE Trans. Signal Process.},
+    author = {Wang, M. and Nehorai, A.},
+    month = feb,
+    year = {2017},
+    pages = {933--946}
+}
+@inproceedings{wang_performance_2017,
+    title = {Performance analysis of coarray-based {MUSIC} and the {Cram\'{e}r}-{Rao} bound},
+    doi = {10.1109/ICASSP.2017.7952719},
+    booktitle = {2017 {IEEE} {International} {Conference} on {Acoustics}, {Speech} and {Signal} {Processing} ({ICASSP})},
+    author = {Wang, M. and Zhang, Z. and Nehorai, A.},
+    month = mar,
+    year = {2017},
+    pages = {3061--3065}
+}
+```

examples/experiments/coarrays_music_crb/sim_efficiency.m

@@ -5,7 +5,7 @@
 
 clear(); close all;
 
-n_snapshots = 1;
+n_snapshots = 1; % normalized
 wavelength = 1;
 d_0 = wavelength / 2;
 power_source = 1;

examples/experiments/coarrays_music_crb/sim_mse_acc_2d.m

@@ -11,7 +11,7 @@
 %  =====
 
 use_simple_rss = true;
-wavelength = 1;
+wavelength = 1; % normalized
 d_0 = wavelength / 2;
 doas = linspace(-pi*6/16, pi*5/16, 11);
 n_doas = length(doas);

examples/experiments/coarrays_music_crb/sim_resolution.m

@@ -10,7 +10,7 @@
 clear(); close all;
 
 rv_mode = 'SS';
-wavelength = 1;
+wavelength = 1; % normalized
 d_0 = wavelength / 2;
 power_source = 1;
 power_noise = 1;

examples/experiments/location_errors/README.md

@@ -0,0 +1,20 @@
+## About
+
+This folder contains scripts for analyzing the performance of SS-MUSIC and DA-MUSIC for various sparse linear arrays in the presence of sensor location errors. You can use these scripts to produce figures similar to those in my following paper:
+
+```
+@article{wang_performance_2018,
+    title = {Performance Analysis of Coarray-Based {MUSIC} in the Presence of Sensor Location Errors},
+    volume = {66},
+    issn = {1053-587X},
+    doi = {10.1109/TSP.2018.2824283},
+    number = {12},
+    journal = {IEEE Transactions on Signal Processing},
+    author = {Wang, M. and Zhang, Z. and Nehorai, A.},
+    month = jun,
+    year = {2018},
+    pages = {3074--3085}
+}
+```
+
+Note that in the paper we use the term "location errors". In the scripts we use "position errors", which has exactly the same meaning.

examples/experiments/location_errors/crb_uc_sto_pos_err.m

@@ -0,0 +1,83 @@
+function [CRB, FIM] = crb_uc_sto_pos_err(design, wavelength, doas, p, noise_var, snapshot_count, pos_err_mask)
+%CRB_UC_STO_POS_ERR Stochastic CRB for joint estimation of DOA parameters and
+%the deterministic location errors with the assumption that the sources are
+%uncorrelated.
+%Inputs:
+%   design - Array design.
+%   wavelength - Wavelength.
+%   doas - DOA vector.
+%   p - Source power vector (diagonals of the source covariance matrix).
+%       If all sources share the same power, you can just pass in a scalar.
+%       You can also pass in a diagonal matrix, whose diagonals will be
+%       extracted automatically.
+%   noise_var - Noise power.
+%   snapshot_count - (Optional) number of snapshots. Default is one.
+%   pos_err_mask - (Optional) A 2xn or 1xn boolean mask specifying the position
+%                  error parameters to be considered. If omitted, the masks will
+%                  be constructed from the non-zero elements from
+%                  design.position_errors.
+%Outputs:
+%   CRB - The CRB matrix.
+%   FIM - The full FIM.
+if design.dim ~= 1
+    error('1D array expected.');
+end
+if nargin <= 5
+    snapshot_count = 1;
+end
+m = design.element_count;
+% Parse position error masks.
+if nargin <= 6
+    % No mask given. Generate from the position errors in the given design.
+    [mask_u, mask_v] = get_pos_err_mask_from_design(design);
+else
+    if size(pos_err_mask, 1) == 1
+        mask_u = pos_err_mask;
+        mask_v = true(1, m);
+    else
+        mask_u = pos_err_mask(1,:);
+        mask_v = pos_err_mask(2,:);
+    end
+end
+% Fall back to the original CRB formula if no position errors are present.
+if ~any(mask_u) && ~any(mask_v)
+    [CRB, FIM] = crb_uc_sto_1d(design, wavelength, doas, p, noise_var, snapshot_count);
+    return;
+end
+% Calculate the CRB according to (33)
+k = length(doas);
+p = unify_source_power_vector(p, k);
+compute_u = any(mask_u);
+compute_v = any(mask_v);
+I = eye(m);
+L1 = I(:, mask_u);
+L2 = I(:, mask_v);
+% Note that we need to use the perturbed steering matrix here.
+[A, DA] = steering_matrix(design, wavelength, doas);
+R = A * bsxfun(@times, p, A') + noise_var * eye(m);
+% Compute each blocks.
+% Note: the Kronecker product may not be efficient for large arrays.
+M_theta = (khatri_rao(conj(DA), A) + khatri_rao(conj(A), DA)) * diag(p);
+M_p = khatri_rao(conj(A), A);
+if compute_u
+    Au = 1j*2*pi/wavelength * A * diag(sin(doas));
+    APAu = A * diag(p) * Au';
+    M_u = khatri_rao(conj(APAu) * L1, L1) + khatri_rao(L1, APAu * L1); 
+else
+    M_u = zeros(m*m, 0);
+end
+if compute_v
+    Av = 1j*2*pi/wavelength * A * diag(cos(doas));
+    APAv = A * diag(p) * Av';
+    M_v = khatri_rao(conj(APAv) * L2, L2) + khatri_rao(L2, APAv * L2); 
+else
+    M_v = zeros(m*m, 0);
+end
+M_n = I(:);
+M = [M_theta M_p M_u M_v M_n];
+FIM = M' / kron(R.', R) * M;
+FIM = 0.5 * real(FIM + FIM') * snapshot_count;
+CRB = eye(size(FIM)) / FIM;
+CRB = CRB(1:k, 1:k);
+end
+

examples/experiments/location_errors/ecov_perturbed_coarray_music_1d.m

@@ -0,0 +1,192 @@
+function C = ecov_perturbed_coarray_music_1d(design, wavelength, doas, p, noise_var, snapshot_count, error_stat, mode)
+%ECOV_PERTURBED_COARRAY_MUSIC_1D Asymptotic covariance matrix of the
+%estimation errors of the coarray-based MUSIC algorithm (SS-MUSIC and
+%DA-MUSIC) in the presence of small model errors.
+%Inputs:
+%   design - Array design.
+%   wavelength - Wavelength.
+%   doas - DOA vector in radians.
+%   p - Source power vector (diagonals of the source covariance matrix).
+%       If all sources share the same power, you can just pass in a scalar.
+%       You can also pass in a diagonal matrix, whose diagonals will be
+%       extracted automatically.
+%   noise_var - Noise power.
+%   snapshot_count - (Optional) number of snapshots. Default is one.
+%   error_stat - A struct describing the statistical properties of the
+%       model errors. The following fields are supported.
+%           'PositionErrorMask' - (Optional) A boolean matrix or vector
+%               specifying which position error parameter is present.
+%               (1) If an 1xM vector is provided, only perturbations along
+%               the x-axis (along the array) will be considered. Set the
+%               m-th element of this vector to true if the m-th element is
+%               perturbed.
+%               (2) If a 2xM vector is provided, perturbations along both
+%               the x-axis and the y-axis (perpendicular to the array) will
+%               be considered. The first row corresponds to the x-axis and
+%               the second row corresponds to the y-axis. For example,
+%               setting the (2,m)-th element of the matrix to true if the
+%               m-th element is perturbed along the y-axis.
+%               (3) If not specified, the following default value will be
+%               used (assuming the first element as the reference element):
+%                   [false true true ... true;...
+%                    false true true ... true]
+%            'PositionErrorCov' - The covariance matrix of the position
+%               errors. The dimension of this matrix must match the number
+%               of true values in PositionErrorMask.
+%   mode - 'Full' or 'DiagonalsOnly'. Default value if 'Full'.
+%Output:
+%   C - If mode is set to 'Full', returns the full asymptotic error
+%       covariance matrix. If mode is set to 'DiagonalsOnly', will only
+%       return a vector of the diagonals. This is useful when you only need
+%       the asymptotic variances for each DOA.
+if design.dim ~= 1
+    error('1D array expected.');
+end
+if nargin <= 6
+    mode = 'full';
+    if nargin <= 5
+        snapshot_count = 1;
+    end
+end
+% Fall back to the perturbation free version if no model errors are present.
+if isempty(error_stat) || isempty(fieldnames(error_stat))
+    C = ecov_coarray_music_1d(design, wavelength, doas, p, noise_var, snapshot_count, mode);
+    return;
+end
+
+m = design.element_count;
+k = length(doas);
+p = unify_source_power_vector(p, k);
+% Generate the selection matrix.
+F = coarray_selection_matrix_1d(design);
+% Check source number.
+m_v = (size(F,1)+1)/2;
+if k >= m_v
+    error('Too many sources.');
+end
+
+% Compute analytical MSE
+% ======================
+
+% Nominal R
+design0 = design;
+design0.mutual_coupling = [];
+design0.position_errors = [];
+design0.gain_errors = [];
+design0.phase_errors = [];
+A0 = steering_matrix(design0, wavelength, doas);
+R0 = A0 * bsxfun(@times, p, A0') + eye(m) * noise_var;
+
+% Coarray
+arr_virtual = ula_1d(m_v, wavelength/2);
+[Av, DAv] = steering_matrix(arr_virtual, wavelength, doas);
+Rss = Av * bsxfun(@times, p, Av') + eye(m_v) * noise_var;
+G = zeros(m_v*m_v, 2*m_v-1);
+
+% Concatenated subarray selection matrix
+for ii = 1:m_v
+    G(((ii-1)*m_v+1):ii*m_v, :) = [zeros(m_v, m_v-ii) eye(m_v) zeros(m_v, ii-1)];
+end
+% Eigendecomposition of the ideal augmented covariance matrix
+[E, ~] = eig(0.5*(Rss + Rss'), 'vector');
+En = E(:,1:end-k);
+% Evalute xi_k / p_k
+Xi_g = zeros(m*m, k);
+gammas = zeros(k, 1);
+pinv_Av = pinv(Av);
+for kk = 1:k
+    d_ak = DAv(:,kk);
+    alpha_k = -pinv_Av(kk,:).';
+    beta_k = En*En'*d_ak;
+    gammas(kk) = real((En'*d_ak)'*(En'*d_ak));
+    Xi_g(:,kk) = F'*G'*kron(beta_k, alpha_k) / gammas(kk) / p(kk);
+end
+% Add model errors' contributions
+% Note: Only sensor location errors are implemented.
+n_perturb = 0; % number of types of perturbations
+B = cell(3,1); % linearized perturbation effect matrix, \Delta r = B\delta
+S = cell(3,1); % covariance matrix of the perturbation parameters
+if isfield(error_stat, 'PositionErrorCov')
+    if isfield(error_stat, 'PositionErrorMask')
+        [nr_mask, nc_mask] = size(error_stat.PositionErrorMask);
+        % dimension check
+        if nc_mask ~= m
+            error('The number of columns of ''PositionErrorMask'' does not match the number of array elements.');
+        end
+        if nr_mask < 1 || nr_mask > 2
+            error('The number of rows of ''PositionErrorMask'' should be either 1 or 2.');
+        end
+        % convert
+        if nr_mask == 1
+            cur_mask = [error_stat.PositionErrorMask false(1, m)];
+        else
+            cur_mask = [error_stat.PositionErrorMask(1,:) error_stat.PositionErrorMask(2,:)];
+        end
+    else
+        % assuming the first element as the reference element
+        cur_mask = true(1,2*m);
+        cur_mask([1 m+1]) = false;
+    end
+    n_perturb = n_perturb + 1;
+    n_pos_err = length(find(cur_mask));
+    cur_S = error_stat.PositionErrorCov;
+    if isscalar(cur_S)
+        % scalar case
+        cur_S = cur_c * eye(n_pos_err);
+    else
+        if any(n_pos_err ~= size(error_stat.PositionErrorCov))
+            error('The dimension %d of ''PositionErrorCov'' does not match the number of position error parameters %d.',...
+                size(error_stat.PositionErrorCov, 1), n_pos_err);
+        end
+    end
+    S{n_perturb} = cur_S;
+
+    % alternative way of constructing B
+%     Ts = A*bsxfun(@times, sin(doas') .* p, A');
+%     Tc = A*bsxfun(@times, cos(doas') .* p, A');
+%     Ms1 = zeros(m*m, m);
+%     Ms2 = zeros(m*m, m);
+%     Mc1 = zeros(m*m, m);
+%     Mc2 = zeros(m*m, m);
+%     for ii = 1:m
+%         Ms1(m*(ii-1)+1:m*ii,ii) = Ts(:,ii);
+%         Ms2(m*(ii-1)+1:m*ii,:) = diag(Ts(:,ii));
+%         Mc1(m*(ii-1)+1:m*ii,ii) = Tc(:,ii);
+%         Mc2(m*(ii-1)+1:m*ii,:) = diag(Tc(:,ii));
+%     end
+%     cur_B = [(Ms1 - Ms2) (Mc1 - Mc2)];
+%     cur_B = (2j * pi / wavelength) * cur_B(:,cur_mask);
+%     B{n_perturb} = cur_B;
+
+    % Using equation (21a) (21b)
+    Au = 1j*2*pi / wavelength * A0 * diag(sin(doas));
+    Av = 1j*2*pi / wavelength * A0 * diag(cos(doas));
+    P = diag(p);
+    Bu = khatri_rao(eye(m), A0 * P * Au') + khatri_rao(conj(A0 * P * Au'), eye(m));
+    Bv = khatri_rao(eye(m), A0 * P * Av') + khatri_rao(conj(A0 * P * Av'), eye(m));
+    cur_B = [Bu Bv];
+    B{n_perturb} = cur_B(:, cur_mask);
+end
+
+% Evaluate cov
+switch lower(mode)
+    case 'diagonalsonly'
+        C = zeros(k,1);
+        for kk = 1:k
+            C(kk) = real(Xi_g(:,kk)' * kron(R0, R0.') * Xi_g(:,kk)) / snapshot_count;
+            for pp = 1:n_perturb
+                T = real(Xi_g(:,kk).' * B{pp});
+                C(kk) = C(kk) + T * S{pp} * T';
+            end
+        end
+    case 'full'
+        C = real(Xi_g' * kron(R0, R0.') * Xi_g / snapshot_count);
+        for pp = 1:n_perturb
+            T = real(Xi_g.' * B{pp});
+            C = C + T * S{pp} * T';
+        end
+    otherwise
+        error('Invalid mode.');
+end
+end
+