diff --git a/.gitignore b/.gitignore
index f8ccc48..a9d86b7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,4 @@
 *.asv
 *.fig
+.vscode
 experimental/
\ No newline at end of file
diff --git a/estimator/mle_sto_uc_1d.m b/estimator/mle_sto_uc_1d.m
index 9ead0f1..896901a 100644
--- a/estimator/mle_sto_uc_1d.m
+++ b/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));
diff --git a/examples/experiments/coarrays_music_crb/README.md b/examples/experiments/coarrays_music_crb/README.md
new file mode 100644
index 0000000..e1e0caf
--- /dev/null
+++ b/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}
+}
+```
diff --git a/examples/experiments/coarrays_music_crb/sim_efficiency.m b/examples/experiments/coarrays_music_crb/sim_efficiency.m
index 38c9307..eea35f8 100644
--- a/examples/experiments/coarrays_music_crb/sim_efficiency.m
+++ b/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;
diff --git a/examples/experiments/coarrays_music_crb/sim_mse_acc_2d.m b/examples/experiments/coarrays_music_crb/sim_mse_acc_2d.m
index e4bf9c1..a0bd638 100644
--- a/examples/experiments/coarrays_music_crb/sim_mse_acc_2d.m
+++ b/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);
diff --git a/examples/experiments/coarrays_music_crb/sim_resolution.m b/examples/experiments/coarrays_music_crb/sim_resolution.m
index 38fba17..1407c5d 100644
--- a/examples/experiments/coarrays_music_crb/sim_resolution.m
+++ b/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;
diff --git a/examples/experiments/location_errors/README.md b/examples/experiments/location_errors/README.md
new file mode 100644
index 0000000..4598cf4
--- /dev/null
+++ b/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.
diff --git a/examples/experiments/location_errors/crb_uc_sto_pos_err.m b/examples/experiments/location_errors/crb_uc_sto_pos_err.m
new file mode 100644
index 0000000..886b548
--- /dev/null
+++ b/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
+
diff --git a/examples/experiments/location_errors/ecov_perturbed_coarray_music_1d.m b/examples/experiments/location_errors/ecov_perturbed_coarray_music_1d.m
new file mode 100644
index 0000000..33cabbf
--- /dev/null
+++ b/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
+
diff --git a/examples/experiments/location_errors/gen_pos_err.m b/examples/experiments/location_errors/gen_pos_err.m
new file mode 100644
index 0000000..399f1e7
--- /dev/null
+++ b/examples/experiments/location_errors/gen_pos_err.m
@@ -0,0 +1,51 @@
+function [err, C] = gen_pos_err(std, n, type, absolute)
+%GEN_POS_ERR Generates position errors. We first generate i.i.d. position errors
+%for each sensor. We use the first sensor as the reference sensor and substract
+%the the position error of the first sensor from the remaining sensors.
+%Therefore, the resulting covariance matrix is not diagonal.
+%Inputs:
+%   std - Standard deviation of the position errors.
+%   n - Number of sensors.
+%   type - 'Gaussian' or 'Uniform'.
+%   absolute - If set to true, will not use the first sensor as the reference
+%              sensor and the position errors will be i.i.d.
+%Outputs:
+%   err - 2xn vector of position errors.
+%   C - Covariance matrix. (2n-2) x (2n-2) if absolute is false, (2n) x (2n) if
+%       absolute is true.
+if nargin < 4
+    absolute = false;
+end
+switch lower(type)
+    case 'gaussian'
+        err = randn(2, n)*std;
+        if nargout > 1
+            if absolute
+                C = std^2*eye(n + n);
+            else
+                C = std^2*(eye(n - 1) + ones(n - 1));
+                C = blkdiag(C, C);
+            end
+        end
+    case 'uniform'
+        b = sqrt(3)*std;
+        err = unifrnd(-b, b, 2, n);
+        if nargout > 1
+            if absolute
+                C = std^2*eye(n + n);
+            else
+                C = std^2*(eye(n - 1) + ones(n - 1));
+                C = blkdiag(C, C);
+            end
+        end
+    otherwise
+        error('Unknow type "%s".', type);
+end
+if ~absolute
+    err(1,2:end) = err(1,2:end) - err(1,1);
+    err(2,2:end) = err(2,2:end) - err(2,1);
+    err(1,1) = 0;
+    err(2,1) = 0;
+end
+end
+
diff --git a/examples/experiments/location_errors/get_design_set.m b/examples/experiments/location_errors/get_design_set.m
new file mode 100644
index 0000000..383d42c
--- /dev/null
+++ b/examples/experiments/location_errors/get_design_set.m
@@ -0,0 +1,40 @@
+function designs = get_design_set(name, d_0)
+%GET_DESIGN_SET Gets a set of arrays designs for simulations.
+%Inputs:
+%   name - 'SameAperture': four designs with the same aperture.
+%          'SameNumSensor': four designs with the same number of sensors.
+%          'Nested10': seven nested arrays with the same number of sensors.
+%   d_0 - Minimal inter-element spacing.
+%Output:
+%   designs - A cell array of array designs.
+switch lower(name)
+    case 'sameaperture'
+        % same aperture
+        designs = { ...
+            design_array_1d('coprime', [2 3], d_0, '2M', 'Co-prime (2,3)') ...
+            design_array_1d('custom', [0 1 2 6 9]*d_0, d_0, 'MRA 5') ...
+            design_array_1d('nested', [1 5], d_0, 'Nested (1,5)') ...
+            design_array_1d('nested', [4 2], d_0, 'Nested (4,2)') ...
+        };
+    case 'samenumsensor'
+        designs = { ...
+            design_array_1d('coprime', [3 5], d_0, '2M', 'Co-prime (3,5)') ...
+            design_array_1d('custom', [0 1 4 10 16 22 28 30 33 35]*d_0, d_0, 'MRA 10') ...
+            design_array_1d('nested', [4 6], d_0, 'Nested (4,6)') ...
+            design_array_1d('nested', [3 7], d_0, 'Nested (3,7)') ...
+        };
+    case 'nested10'
+        designs = { ...
+            design_array_1d('nested', [2 8], d_0, 'Nested (2,8)') ...
+            design_array_1d('nested', [3 7], d_0, 'Nested (3,7)') ...
+            design_array_1d('nested', [4 6], d_0, 'Nested (4,6)') ...
+            design_array_1d('nested', [5 5], d_0, 'Nested (5,5)') ...
+            design_array_1d('nested', [6 4], d_0, 'Nested (6,4)') ...
+            design_array_1d('nested', [7 3], d_0, 'Nested (7,3)') ...
+            design_array_1d('nested', [8 2], d_0, 'Nested (8,2)') ...
+        };
+    otherwise
+        error('Unknown design set ''%s''.', name);
+end
+end
+
diff --git a/examples/experiments/location_errors/get_pos_err_mask_from_design.m b/examples/experiments/location_errors/get_pos_err_mask_from_design.m
new file mode 100644
index 0000000..7b86584
--- /dev/null
+++ b/examples/experiments/location_errors/get_pos_err_mask_from_design.m
@@ -0,0 +1,24 @@
+function [mask_u, mask_v] = get_pos_err_mask_from_design(design)
+%GET_POS_ERR_MASK_FROM_DESIGN Gets mask vectors from the position_errors field
+%in the given design. The masks are generated by identifying non-zero elements
+%in design.position_errors.
+%Input:
+%   design - Array design.
+%Outputs:
+%   mask_u - Mask vector for the position errors along the x-axis.
+%   mask_v - Mask vector for the position errors along the y-axis.
+
+if ~isfield(design, 'position_errors') || isempty(design.position_errors)
+    mask_u = false(1, design.element_count);
+    mask_v = false(1, design.element_count);
+    return
+end
+
+mask_u = design.position_errors(1,:) ~= 0;
+if size(design.position_errors, 1) == 1
+    mask_v = false(1, design.element_count);
+else
+    mask_v = design.position_errors(2,:) ~= 0;
+end
+
+end
\ No newline at end of file
diff --git a/examples/experiments/location_errors/mle_uc_sto_pos_err.m b/examples/experiments/location_errors/mle_uc_sto_pos_err.m
new file mode 100644
index 0000000..03f93b7
--- /dev/null
+++ b/examples/experiments/location_errors/mle_uc_sto_pos_err.m
@@ -0,0 +1,107 @@
+function [doas, p, u, v, noise_var] = mle_uc_sto_pos_err(R, n_doas, design, wavelength, x0, pos_err_mask)
+%MLE_UC_STO_POS_ERR Stochastic maximum-likelihood estimator with the
+%uncorrelated source assumption for jointly estimate the DOA parameters and
+%the sensor location errors.
+%Inputs:
+%   R - Sample covariance matrix.
+%   n_doas - Number of sources.
+%   design - Array design.
+%   wavelength - Wavelength.
+%   x_0 - (Optional) an vector storing the initial guess in the following order:
+%         [doas; p; u; v; noise_var]. If absent, an initial guess will be
+%         obtained using SS-MUSIC. The input DOAs must be in radians.
+%   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:
+%   doas - DOA estimates in radians.
+%   p - Estimates of the source powers.
+%   u - Estimates of the position errors along the x-axis.
+%   v - Estimates of the position errors along the y-axis.
+%   noise_var - Estimate of the noise variance.
+if nargin < 5
+    x0 = [];
+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
+% Find number of unknown location errors.
+n_u = length(find(mask_u));
+n_v = length(find(mask_v));
+n_unknown = n_doas + n_doas + n_u + n_v + 1;
+if ~isempty(x0) && length(x0) ~= n_unknown
+    error('The vector of the initial guess has incorrect length.');
+end
+design_nominal = design;
+design_nominal.position_errors = [];
+% Compute the initial guess.
+if isempty(x0)
+    x0 = zeros(n_unknown, 1);
+    % Use SS-MUSIC to obtain the initial guess of the DOAs.
+    [Rss, ~] = virtual_ula_cov_1d(design, R, 'SS');
+    sp = rmusic_1d(Rss, n_doas, 2 * pi * design.element_spacing / wavelength);
+    x0(1:n_doas) = sp.x_est;
+    % Use least squares to obtain the initial guess of the source powers
+    % and noise power.
+    A_est = steering_matrix(design_nominal, wavelength, sp.x_est);
+    pn_est = real(linsolve([khatri_rao(conj(A_est), A_est) reshape(eye(m), [], 1)], R(:)));
+    x0(n_doas+1:n_doas+n_doas) = pn_est(1:end-1);
+    x0(end) = pn_est(end);
+    % Position errors remain zeros.
+end
+% NLL wrapper - decompose x into difference components and pass them to the NLL.
+f = @(x) nll(R, design_nominal, wavelength, ...
+    x(1:n_doas), x(n_doas+1:n_doas+n_doas), ...
+    x(n_doas+n_doas+1:n_doas+n_doas+n_u), ...
+    x(n_doas+n_doas+n_u+1:n_doas+n_doas+n_u+n_v), ...
+    mask_u, mask_v, x(end));
+% Constraints.
+lb = [ones(n_doas, 1) * (-pi/2);...
+      zeros(n_doas, 1);...
+      ones(n_u + n_v, 1) * (-wavelength);...
+      0];
+ub = [ones(n_doas, 1) * pi/2;...
+      inf(n_doas, 1);...
+      ones(n_u + n_v, 1) * wavelength;...
+      inf];
+% Call fmincon. This will take a long time.
+options = optimoptions('fmincon', 'MaxFunctionEvaluations', 50000, 'Display', 'none');
+sol = fmincon(f, x0, [], [], [], [], lb, ub, [], options);
+% Output the solution.
+doas = sol(1:n_doas);
+p = sol(n_doas+1:n_doas+n_doas);
+u = sol(n_doas+n_doas+1:n_doas+n_doas+n_u);
+v = sol(n_doas+n_doas+n_u+1:n_doas+n_doas+n_u+n_v);
+noise_var = sol(end);
+end
+
+function obj = nll(R, design, wavelength, doas, p, u, v, mask_u, mask_v, noise_var)
+%NLL Negative log-likelihood function to be optimized.
+pos_errs = zeros(2, design.element_count);
+pos_errs(1,mask_u) = u;
+pos_errs(2,mask_v) = v;
+design_perturbed = design;
+design_perturbed.position_errors = pos_errs;
+A = steering_matrix(design_perturbed, wavelength, doas);
+S = A * diag(p) * A' + noise_var*eye(design.element_count);
+S = 0.5 * (S + S');
+[L,p] = chol(S);
+if p > 0
+    obj = inf;
+    return;
+end
+logdetS = sum(log(diag(L)));
+obj = logdetS + sum(real(diag(R / S)));
+end
diff --git a/examples/experiments/location_errors/sim_crb_vs_mle_fixed_pos_err.m b/examples/experiments/location_errors/sim_crb_vs_mle_fixed_pos_err.m
new file mode 100644
index 0000000..0cb5594
--- /dev/null
+++ b/examples/experiments/location_errors/sim_crb_vs_mle_fixed_pos_err.m
@@ -0,0 +1,85 @@
+% Demonstrates the CRB obtained from (33) is indeed attainable by the stochastic
+% MLE.
+% Produces figures similar to Fig. 9.
+
+%% Configure
+clear();
+rng(42);
+
+wavelength = 1; % normalized
+d_0 = wavelength / 2;
+design_set = 'SameNumSensor'; % must have the same number of sensors!
+designs = get_design_set(design_set, d_0);
+n_designs = length(designs);
+n_snapshots = 5000;
+n_doas = 11;
+% Note: here these numbers are reduces because the MLE takes a long time to
+% compute. Consider increasing them to obtain a more stable result.
+n_repeats = 10;
+n_params = 10;
+doas = linspace(-pi/3, pi/3, n_doas);
+
+power_source = 1;
+params_SNR = linspace(-20, 20, n_params);
+
+% generate fixed position error vectors
+pos_err = randn(2, designs{1}.element_count) * (0.1 * d_0);
+pos_err(:,[1 2]) = 0;
+mask_u = pos_err(1,:) ~= 0;
+mask_v = pos_err(2,:) ~= 0;
+pos_err_mask = [mask_u; mask_v];
+
+%% Run
+mses = zeros(n_designs, n_params, n_repeats);
+crbs = zeros(n_designs, n_params);
+crbs_nominal = zeros(n_designs, n_params);
+
+for dd = 1:n_designs
+    design_nominal = designs{dd};
+    design = design_nominal;
+    design.position_errors = pos_err(:,1:design.element_count);
+    fprintf('Running simulation for ''%s''\n', design.name);
+    
+    % Remember to remove the progressbar calls when using parfor.
+    progressbar('reset', n_params * n_repeats);
+    for ii = 1:n_params
+        SNR = params_SNR(ii);
+        power_noise = power_source * 10^(-SNR/10);
+        for rr = 1:n_repeats
+            [~, R] = snapshot_gen_sto(design, doas, wavelength, n_snapshots, power_noise, power_source);
+            [doa_est, p_est, u_est, v_est, noise_est] = mle_uc_sto_pos_err(R, n_doas, design, wavelength, [], pos_err_mask);
+            mses(dd, ii, rr) = mean((doa_est(:) - doas(:)).^2);
+            progressbar('advance');
+        end
+        CRB = crb_uc_sto_pos_err(design, wavelength, doas, power_source, power_noise, n_snapshots, pos_err_mask);
+        crbs(dd, ii) = mean(diag(CRB));
+        CRB_nominal = crb_uc_sto_1d(design_nominal, wavelength, doas, power_source, power_noise, n_snapshots);
+        crbs_nominal(dd, ii) = mean(diag(CRB_nominal));
+    end
+    progressbar('end');
+end
+    
+%% Plot
+% The MLE obtained from fmincon may produce very large errors, leading to
+% outliers which will greatly affect the empirical MSE. Set the following
+% variable to true to suppress DOA estimates that are too large.
+suppress_outliers = false;
+processed_mses = mses;
+if suppress_outliers
+    max_mse = ((max(doas) - min(doas)) / (n_doas - 1) / 2)^2;
+    processed_mses(processed_mses > max_mse) = nan;
+end
+mses_final = nanmean(processed_mses, 3);
+
+for dd = 1:n_designs
+    hf = figure;
+    semilogy(params_SNR, mses_final(dd,:), 'kx', ...
+        params_SNR, crbs(dd,:), 'k--', params_SNR, crbs_nominal(dd,:), 'k-.');
+    legend('MSE', 'CRB', 'CRB (location-error free)');
+    axis([-inf inf -inf 1]);
+    xlabel('SNR (dB)'); ylabel('MSE (rad^2)');
+    title(designs{dd}.name);
+    fig_pos = get(hf, 'Position');
+    set(hf, 'Position', [fig_pos(1) fig_pos(2) 400 300]);
+    set(hf, 'Color', 'White');
+end
diff --git a/examples/experiments/location_errors/sim_error_analysis_det_full.m b/examples/experiments/location_errors/sim_error_analysis_det_full.m
new file mode 100644
index 0000000..085821c
--- /dev/null
+++ b/examples/experiments/location_errors/sim_error_analysis_det_full.m
@@ -0,0 +1,83 @@
+% Deterministic error model.
+% Full verification.
+% Produces figures similar to Fig. 3.
+
+%% Configure
+clear();
+
+wavelength = 1; % normalized
+d_0 = wavelength / 2;
+designs = get_design_set('SameNumSensor', d_0);
+
+n_designs = length(designs);
+n_doas = 11;
+doas = linspace(-pi/3, pi/3, n_doas);
+
+source_power = 1;
+noise_power = 1;
+
+perturb_type = 'gaussian';
+n_params = 20; % 20 is used in the paper.
+max_perturb = 0.10*d_0;
+max_n_snapshots = 1000;
+params_perturb = linspace(0, max_perturb, n_params);
+params_n_snapshots = round(linspace(100, max_n_snapshots, n_params));
+
+n_repeats = 1000;
+%% Run
+mse_em = zeros(n_designs, n_params, n_params, n_repeats);
+mse_an = zeros(n_designs, n_params, n_params);
+
+fprintf('Varying position error std and number of snapshots:\n');
+fprintf('    DOAs: [%s]\n', num2str(doas, '%f '));
+fprintf('    SNR: %.1f\n', 10*log10(source_power / noise_power));
+fprintf('    Repeat: %d\n', n_repeats);
+
+for dd = 1:n_designs
+    cur_design = designs{dd};
+    fprintf('Running simulations for %s:\n', cur_design.name);
+    progressbar('reset', n_params^2);
+    for ii = 1:n_params
+        n_snapshots = params_n_snapshots(ii);
+        for kk = 1:n_params
+            perturb_std = params_perturb(kk);
+            % collect empirical results
+            for rr = 1:n_repeats
+                pos_err = gen_pos_err(perturb_std, cur_design.element_count, perturb_type);
+                perturbed_design = cur_design;
+                perturbed_design.position_errors = pos_err;
+                [~, R] = snapshot_gen_sto(perturbed_design, doas, wavelength,...
+                            n_snapshots, noise_power, source_power);
+                [Rv, ~, ~] = virtual_ula_cov_1d(cur_design, R, 'SS');
+                sp = rmusic_1d(Rv, n_doas, 2*pi * cur_design.element_spacing / wavelength);
+                mse_em(dd,ii,kk,rr) = sum((sp.x_est - doas).^2) / n_doas;
+            end
+            % compute analytical results
+            [~, perturb_cov] = gen_pos_err(perturb_std, cur_design.element_count, perturb_type);
+            error_stat = struct;
+            error_stat.PositionErrorCov = perturb_cov;
+            mse_an(dd,ii,kk) = sum(ecov_perturbed_coarray_music_1d(cur_design, wavelength, ...
+                doas, source_power, noise_power, n_snapshots, error_stat, 'DiagonalsOnly')) / n_doas;
+            progressbar('advance');
+        end
+    end
+    progressbar('end');
+end
+
+%% Plot
+rmse_em = rad2deg(sqrt(mean(mse_em, 4)));
+rmse_an = rad2deg(sqrt(mse_an));
+rel_err = abs(rmse_em - rmse_an) ./ rmse_an;
+
+for dd = 1:n_designs
+    figure;
+    imagesc(params_perturb/d_0, params_n_snapshots, squeeze(rel_err(dd,:,:))); axis ij;
+    xlabel('\delta/d_0');
+    ylabel('Number of Snapshots');
+    title(designs{dd}.name);
+    colormap gray;
+    colormap(flipud(colormap));
+    colorbar;
+    set(gca, 'Clim', [0 0.5]);
+    set(gca,'YDir','normal');
+end
diff --git a/examples/experiments/location_errors/sim_error_analysis_det_perturb.m b/examples/experiments/location_errors/sim_error_analysis_det_perturb.m
new file mode 100644
index 0000000..107e9c0
--- /dev/null
+++ b/examples/experiments/location_errors/sim_error_analysis_det_perturb.m
@@ -0,0 +1,87 @@
+% Deterministic error model.
+% Demonstrates how the MSE various as the perturbation level increases.
+% Produces figures similar to Fig. 4, 5.
+
+%% Configure
+clear();
+
+wavelength = 1; % normalized
+d_0 = wavelength / 2;
+design_set = 'SameAperture'; % can also be 'SameAperture'
+designs = get_design_set(design_set, d_0);
+
+n_designs = length(designs);
+if strcmpi(design_set, 'SameNumSensor')
+    n_doas = 11;
+else
+    n_doas = 6;
+end
+doas = linspace(-pi/3, pi/3, n_doas);
+
+n_snapshots = 1000;
+source_power = 1;
+noise_power = 1;
+
+perturb_type = 'gaussian';
+n_params_perturb = 20;
+max_perturb = 0.08*d_0;
+params_perturb = linspace(0, max_perturb, n_params_perturb);
+n_repeats = 1000;
+
+%% Run
+mse_em = zeros(n_designs, n_params_perturb, n_repeats);
+mse_an = zeros(n_designs, n_params_perturb);
+
+fprintf('Fixed snapshot count, varying pos err std parameters:\n');
+fprintf('    DOAs: [%s]\n', num2str(doas, '%f '));
+fprintf('    Number of snapshots: %d\n', n_snapshots);
+fprintf('    SNR: %.1f\n', 10*log10(source_power / noise_power));
+fprintf('    Perturbation range: [0 %f]d_0\n', max_perturb/d_0);
+fprintf('    Repeat: %d\n', n_repeats);
+
+for dd = 1:n_designs
+    cur_design = designs{dd};
+    fprintf('Running simulations for %s:\n', cur_design.name);
+    progressbar('reset', n_params_perturb);
+    for kk = 1:n_params_perturb
+        % collect empirical results
+        cur_mses = zeros(n_repeats, 1);
+        perturb_std = params_perturb(kk);
+        for rr = 1:n_repeats
+            pos_err = gen_pos_err(perturb_std, cur_design.element_count, perturb_type);
+            perturbed_design = cur_design;
+            perturbed_design.position_errors = pos_err;
+            [~, R] = snapshot_gen_sto(perturbed_design, doas, wavelength, n_snapshots, noise_power, source_power);
+            [Rv, ~, ~] = virtual_ula_cov_1d(cur_design, R, 'SS');
+            sp = rmusic_1d(Rv, n_doas, 2*pi * cur_design.element_spacing / wavelength);
+            mse_em(dd,kk,rr) = sum((sp.x_est - doas).^2) / n_doas;
+        end
+        % compute analytical results
+        [~, perturb_cov] = gen_pos_err(params_perturb(kk), cur_design.element_count, perturb_type);
+        error_stat = struct;
+        error_stat.PositionErrorCov = perturb_cov;
+        mse_an(dd, kk) = sum(ecov_perturbed_coarray_music_1d(cur_design, wavelength, ...
+            doas, source_power, noise_power, n_snapshots, error_stat, 'DiagonalsOnly')) / n_doas;
+        progressbar('advance');
+    end
+    progressbar('end');
+end
+
+%% Plot
+markers = {'x', 'o', 's', 'd', '*', '^', '>', '<', 'h'};
+figure;
+for dd = 1:n_designs
+    hp = plot(params_perturb/d_0, rad2deg(sqrt(mean(squeeze(mse_em(dd,:,:)), 2))), ...
+        ['' markers{mod(dd, length(markers))+1}], ...
+        'DisplayName', [designs{dd}.name ' empirical']);
+    hold on;
+    plot(params_perturb/d_0, rad2deg(sqrt(mse_an(dd,:))), '--', 'Color', get(hp, 'Color'), ...
+        'DisplayName', [designs{dd}.name ' analytical']);
+end
+hold off;
+xlabel('\delta_p/d_0');
+ylabel('RMSE/deg');
+grid on;
+legend('show');
+legend('Location', 'northwest');
+title('RMSE vs Perturbation Level - Deterministic Error Model');
diff --git a/examples/experiments/location_errors/sim_error_analysis_det_snr.m b/examples/experiments/location_errors/sim_error_analysis_det_snr.m
new file mode 100644
index 0000000..b43d6cd
--- /dev/null
+++ b/examples/experiments/location_errors/sim_error_analysis_det_snr.m
@@ -0,0 +1,82 @@
+% Deterministic error model.
+% MSE vs. SNR under different levels of perturbations.
+% Produces figures similar to Fig. 6.
+
+%% Configure
+clear();
+
+wavelength = 1; % normalized
+d_0 = wavelength / 2;
+design = design_array_1d('coprime', [2 3], d_0, '2M', 'Co-prime (2,3)');
+
+n_doas = 6;
+doas = linspace(-pi/3, pi/3, n_doas);
+
+n_snapshots = 5000;
+source_power = 1;
+noise_power = 1;
+
+perturb_type = 'gaussian';
+params_perturb = [0 0.01 0.02 0.04] * d_0;
+n_params_perturb = length(params_perturb);
+max_perturb = max(params_perturb);
+n_params_snr = 20;
+params_snr = linspace(-15, 15, n_params_snr);
+n_repeat = 1000;
+
+%% Run
+
+mse_em_ss = zeros(n_params_perturb, n_params_snr, n_repeat);
+mse_an = zeros(n_params_perturb, n_params_snr);
+
+fprintf('Fixed snapshot count, varying SNR parameters:\n');
+fprintf('    DOAs: [%s]\n', num2str(doas, '%f '));
+fprintf('    Number of snapshots: %d\n', n_snapshots);
+fprintf('    Perturbation range: [0 %f]d_0\n', max_perturb/d_0);
+fprintf('    Repeat: %d\n', n_repeat);
+
+progressbar('reset', n_params_perturb * n_params_snr);
+for kk = 1:n_params_perturb
+    perturb_std = params_perturb(kk);
+    for ss = 1:n_params_snr
+        noise_power = 10^(-params_snr(ss) / 10);
+        % collect empirical results
+        for rr = 1:n_repeat
+            pos_err = gen_pos_err(perturb_std, design.element_count, perturb_type);
+            perturbed_design = design;
+            perturbed_design.position_errors = pos_err;
+            [~, R] = snapshot_gen_sto(perturbed_design, doas, wavelength, n_snapshots, noise_power, source_power);
+            [Rv, ~, ~] = virtual_ula_cov_1d(design, R, 'SS');
+            sp = rmusic_1d(Rv, n_doas, 2*pi * design.element_spacing / wavelength);
+            mse_em_ss(kk, ss, rr) = sum((sp.x_est - doas).^2) / n_doas;
+        end
+        % compute analytical results
+        [~, perturb_cov] = gen_pos_err(perturb_std, design.element_count, perturb_type);
+        error_stat = struct;
+        error_stat.PositionErrorCov = perturb_cov;
+        mse_an(kk, ss) = sum(ecov_perturbed_coarray_music_1d(design, wavelength, ...
+            doas, source_power, noise_power, n_snapshots, error_stat, 'DiagonalsOnly')) / n_doas;
+        progressbar('advance');
+    end
+end
+progressbar('end');
+
+%% Plot
+markers = {'x', 'o', 's', 'd', '*', '^'};
+figure;
+for pp = 1:n_params_perturb
+    perturb_str = sprintf('\\sigma_p/d_0 = %.2f', params_perturb(pp)/d_0);
+    hp = semilogy(params_snr, rad2deg(sqrt(nanmean(squeeze(mse_em_ss(pp,:,:)), 2))), ...
+        ['' markers{mod(pp, length(markers))+1}], ...
+        'DisplayName', ['SSM empirical:' perturb_str]);
+    hold on;
+    semilogy(params_snr, rad2deg(sqrt(mse_an(pp,:))), '--',...
+        'Color', get(hp, 'Color'), ...
+        'DisplayName', ['SSM analytical:' perturb_str]);
+    hold on;
+end
+xlabel('SNR/dB');
+ylabel('RMSE/deg');
+hold off;
+legend('show');
+title('RMSE vs SNR - Deterministic Error Model');
diff --git a/examples/experiments/location_errors/sim_error_analysis_sto.m b/examples/experiments/location_errors/sim_error_analysis_sto.m
new file mode 100644
index 0000000..d8ab7dd
--- /dev/null
+++ b/examples/experiments/location_errors/sim_error_analysis_sto.m
@@ -0,0 +1,82 @@
+% Stochastic error model.
+% RMSE vs number of snapshots.
+% Produces figures similar to Fig. 7 and 8.
+
+%% Configure
+clear(); %close all;
+
+wavelength = 1; % normalized
+d_0 = wavelength / 2;
+design_set = 'SameNumSensor';  % can also be 'SameAperture'
+designs = get_design_set(design_set, d_0);
+n_designs = length(designs);
+if strcmpi(design_set, 'SameNumSensor')
+    n_doas = 11;
+else
+    n_doas = 6;
+end
+doas = linspace(-pi/3, pi/3, n_doas);
+
+source_power = 1;
+noise_power = 1;
+
+perturb_type = 'gaussian';
+perturb_std = 0.1*d_0;
+n_params_n_snapshots = 20;
+params_n_snapshots = floor(linspace(100, 5000, n_params_n_snapshots));
+n_repeats = 50; % set to 5000 in the paper
+
+%% Run
+mse_em = zeros(n_designs, n_params_n_snapshots, n_repeats);
+mse_an = zeros(n_designs, n_params_n_snapshots);
+
+fprintf('Fixed perturbation level, varying number of snapshots:\n');
+fprintf('    DOAs: [%s]\n', num2str(doas, '%f '));
+fprintf('    SNR: %.1f\n', 10*log10(source_power / noise_power));
+fprintf('    # of snapshots: [%d %d]\n', min(params_n_snapshots), max(params_n_snapshots));
+fprintf('    Repeat: %d\n', n_repeats);
+
+for dd = 1:n_designs
+    cur_design = designs{dd};
+    pos_err_info = struct('std', perturb_std, 'mask', ...
+        true(cur_design.element_count, 1), 'reference', 'absolute');
+    fprintf('Running simulations for %s:\n', cur_design.name);
+    progressbar('reset', n_params_n_snapshots);
+    for kk = 1:n_params_n_snapshots
+        n_snapshots = params_n_snapshots(kk);
+        % collect empirical results
+        for rr = 1:n_repeats
+            [~, R] = snapshot_gen_sto_pos_err(cur_design, doas, wavelength, ...
+                n_snapshots, noise_power, source_power, pos_err_info);
+            [Rv, ~, ~] = virtual_ula_cov_1d(cur_design, R, 'SS');
+            sp = rmusic_1d(Rv, n_doas, 2*pi * cur_design.element_spacing / wavelength);
+            mse_em(dd,kk,rr) = sum((sp.x_est - doas).^2) / n_doas;
+        end
+        % compute analytical results
+        c1 = exp(-4 * pi^2 / wavelength^2 * perturb_std^2);
+        new_noise_power = (noise_power + (1 - c1) * n_doas * source_power) / c1;
+        mse_an(dd,kk) = mean(ecov_coarray_music_1d(cur_design, wavelength, doas, ...
+            source_power, new_noise_power, n_snapshots, 'DiagonalsOnly'));
+        progressbar('advance');
+    end
+    progressbar('end');
+end
+
+%% Plot
+markers = {'x', 'o', 's', 'd', '*', '^'};
+hf = figure;
+for dd = 1:n_designs
+    hp = semilogy(params_n_snapshots, rad2deg(sqrt(mean(squeeze(mse_em(dd,:,:)), 2))), ...
+        markers{mod(dd, length(markers))+1}, ...
+        'DisplayName', [designs{dd}.name ' empirical']);
+    hold on;
+    semilogy(params_n_snapshots, rad2deg(sqrt(mse_an(dd,:))), '--',...
+        'Color', get(hp, 'Color'), ...
+        'DisplayName', [designs{dd}.name ' analytical']);
+end
+hold off;
+grid on;
+xlabel('Number of snapshots');
+ylabel('RMSE/deg');
+legend('show', 'Location', 'northeast');
+title('RMSE vs # of Snapshots - Stochastic Error Model');
diff --git a/examples/experiments/location_errors/snapshot_gen_sto_pos_err.m b/examples/experiments/location_errors/snapshot_gen_sto_pos_err.m
new file mode 100644
index 0000000..925468f
--- /dev/null
+++ b/examples/experiments/location_errors/snapshot_gen_sto_pos_err.m
@@ -0,0 +1,79 @@
+function [X, R, S] = snapshot_gen_sto_pos_err(design, doas, wavelength, t, ncov, scov, pos_err_info)
+%SNAPSHOT_GEN_STO_POS_ERR Generates snapshots for the stochastic model
+%with 2D sensor location errors changes every snapshot.
+%Syntax:
+%   X = SNAPSHOT_GEN_STO_POS_ERR_2D(design, doas, wavelength[, t, ncov, scov, pos_err_info]);
+%   [X, R] = SNAPSHOT_GEN_STO_POS_ERR_2D(design, doas, wavelength[, t, ncov, scov, pos_err_info]);
+%   [X, R, S] = SNAPSHOT_GEN_STO_POS_ERR_2D(design, doas, wavelength[, t, ncov, scov, pos_err_info]);
+%Inputs:
+%   design - Array design.
+%   doas - DOA vector. For 2D DOAs, each column represents a DOA pair.
+%   wavelength - Wavelength.
+%   t - Number of snapshots.
+%   ncov - Covariance matrix of the additive complex circular-symmetric
+%          Gaussian noise. Can be a scalar, vector (for uncorrelated noise
+%          with different powers), or a matrix.
+%   scov - Covariance matrix of the source signals. Can be a scalar, vector
+%          (for uncorrelated sources with different powers), or a matrix.
+%   pos_err_info - A struct containing the information about position
+%                  errors:
+%                   'std' - Standard deviation.
+%                   'mask' - Specifies which sensor has location errors.
+%                   'reference' - 'Absolute' or 'First'.
+%Outputs:
+%   X - Snapshots, where each columns is a single snapshot.
+%   R - Sample covariance matrix (averaged by the number of snapshots).
+%   S - A source_count x snapshot_count matrix consists of source signal
+%       vectors.
+if nargin <= 6
+    pos_err_info = struct('std', 0, 'mask', true(design.element_count, 1), ...
+        'reference', 'Absolute');
+end
+if nargin <= 5
+    scov = 1;
+end
+if nargin <= 4
+    ncov = 1;
+end
+if nargin <= 3
+    t = 1;
+end
+m = design.element_count;
+k = length(doas);
+S_internal = gen_ccsg(k, t, scov);
+N = gen_ccsg(m, t, ncov);
+X = zeros(m, t);
+ind = find(pos_err_info.mask);
+n_err = length(ind);
+% lazy comparison here
+ref_first = strcmpi(pos_err_info.reference, 'first');
+for tt = 1:t
+    % generate position errors
+    pos_err = zeros(2, m);
+    pos_err(:,ind) = pos_err_info.std * randn(2, n_err);
+    if ref_first
+        pos_err = pos_err - pos_err(:,1);
+    end
+    design.position_errors = pos_err;
+    A = steering_matrix(design, wavelength, doas);
+    X(:,tt) = A * S_internal(:,tt) + N(:,tt);
+end
+if nargout >= 2
+    R = (X*X')/t;
+    if nargout == 3
+        S = S_internal;
+    end
+end
+end
+
+function X = gen_ccsg(m, n, cov)
+X0 = randn(m, n) + 1j * randn(m, n);
+if isscalar(cov)
+    X = sqrt(cov/2) * X0;
+elseif isvector(cov)
+    X = bsxfun(@times, X0, cov(:));
+else
+    C = sqrtm(cov/2);
+    X = C*X0;
+end
+end
\ No newline at end of file