Array Self-Calibration and Direction Finding in the Presence of Mutual Coupling and Calibration Uncertainties

Since R2026a

Introduction

This example demonstrates an eigenstructure-based method for direction finding in the presence of mutual coupling and calibration uncertainties, using the approach proposed by Friedlander and Weiss [1].

The method enhances the standard MUSIC algorithm by iteratively refining estimates of:

• Direction of arrival (DOA)

• Sensor calibration errors

• Mutual coupling effects between array elements

The technique does not rely on knowledge of source directions and can operate with signals of opportunity.

It applies to both Uniform Linear Arrays (ULAs) and Uniform Circular Arrays (UCAs). These array geometries allow us to simplify the mutual coupling matrix (MCM) based on the array's structure. Specifically, MCMs are assumed to have banded Toeplitz structures in ULAs, and circulant structures in UCAs, enabling us to estimate them directly from observed data.

Background

Subspace-based direction finding algorithms, such as MUSIC (ADD LINK), estimate the DOAs of narrowband sources by finding array steering vectors that are orthogonal to the noise subspace.

However, the presence of mutual coupling and sensor calibration errors distorts the ideal steering vectors, leading to degraded DOA estimates. By including models of the MCM and phase and amplitude calibration errors in the steering vector models, it is possible to improve on subspace-based direction finding algorithms.

Mutual Coupling

Mutual coupling is a phenomenon that occurs in antenna arrays, where a voltage induced in one antenna element generates voltages in nearby elements. This interaction distorts the received signals from their ideal behavior and must be accounted for in array signal processing algorithms.

In this example, mutual coupling is modeled using a mutual coupling matrix (MCM). The MCM is an MxM matrix where M is the number of elements in the antenna array. Each element of this complex-valued matrix represents the amplitude and phase of the signal induced in one antenna element by another.

The diagram below illustrates the interpretation of entries in a single column of the MCM for an array with arbitrary geometry, where $c_{\text{ij}}$ represents coupling from element i to element j:

The MCM would have the following structure:

C = $\left[\begin{array}{ccccc}{c}_{11}& c_{21}& c_{31}& c_{41}& c_{51}\\ c_{12}& c_{22}& c_{32}& c_{42}& c_{52}\\ c_{13}& c_{23}& c_{33}& c_{43}& c_{53}\\ c_{14}& c_{24}& c_{34}& c_{44}& c_{54}\\ c_{15}& c_{25}& c_{35}& c_{45}& c_{55}\end{array}\right]$

To reduce the number of unknowns in the MCM, several assumptions are often made. A common one is symmetry, meaning the mutual coupling from element i to element j is equal to the coupling from j to i. Additional structural assumptions, specific to ULAs and UCAs, are discussed later in this example. In all cases here, we assume the diagonal elements of the MCM are equal to 1, representing each element’s direct response without coupling.

Calibration Errors

In the context of this example, calibration errors refer to phase and amplitude errors that occur within each channel of the array receiver. These errors are present in all phased array systems due to path length differences and variability in the components within each channel. Like mutual coupling, it is important to account for these errors in many signal processing algorithms.

Calibration errors are also represented as an MxM matrix, but the only populated values are along the diagonal, because these calibration errors only impact a single element channel. The calibration matrix for the arbitrary geometry array shown above therefore has the following structure, where each value on the diagonal is a complex number and all other values are zero:

$$G=\left[\begin{array}{ccccc}{g}_{11}& 0& 0& 0& 0\\ 0& g_{22}& 0& 0& 0\\ 0& 0& g_{33}& 0& 0\\ 0& 0& 0& g_{44}& 0\\ 0& 0& 0& 0& g_{55}\end{array}\right]$$

Signal Model

We now present the signal model used in this example. This signal model includes both the mutual coupling and calibration matrices, as well as the array of steering vectors due to the array geometry:

$$X=C*G*A\left(\theta \right)*S+N$$

Where:

• $X$ is the received signal.

• $C$ is the MCM.

• $G$ is the diagonal calibration matrix.

• $A\left(\theta \right)$ is the steering matrix for each signal arriving at the array at angle $\theta$. This matrix contains N columns for N signal sources.

• $S$ is the signal arriving at the array.

• $N$ is the receiver noise.

The typical goal of DOA algorithms is to find$A\left(\theta \right)$. In this example, we demonstrate how to find $A\left(\theta \right)$ in the presence of $C$ and $G$ while improving the estimates of $C$ and $G$.

Traditionally, mutual coupling and calibration matrices are estimated prior to deployment using a variety of techniques. However, these parameters may drift over time due to temperature changes, aging hardware, or reconfiguration of the system. In such cases, online self-calibration becomes necessary to maintain high-fidelity direction finding.

Cost Function

This algorithm builds upon the principles of the MUSIC algorithm [2], but extends it by incorporating the MCM and the calibration matrix into the cost function being minimized. The modified cost function is:

$$J=\sum _{n=1}^N\Vert {\underset{}{\overset{ˆ}{U}}}^H*C*G*a\left({\theta }_n\right)\Vert$$

Where:

• $J$ is the cost function.

• $\underset{}{\overset{ˆ}{U}}$ is an estimate of the noise eigenvectors $U.$

• $a\left({\theta }_n\right)$ are the steering vector for each signal N. These are the columns of $A\left(\theta \right)$.

The core idea is that the product $C*G*a\left({\theta }_n\right)$ should lie in the signal subspace and therefore be orthogonal to the noise subspace represented by $U$, thereby minimizing $J$. This mirrors the philosophy behind MUSIC, but differs in that the cost function here explicitly includes the coupling and calibration effects.

Before estimating our unknown parameters, we must construct $\underset{}{\overset{ˆ}{U}}$. To do so, we collect received signal$X$and estimate the covariance matrix $R_X$. We then perform eigendecomposition on this matrix. The resultant matrix contains N eigenvectors corresponding to the signal subspace, and M - N eigenvectors corresponding to the noise subspace. Assuming that we know N, or that we can estimate it based on the magnitude of the eigenvalues, we can construct a matrix $\underset{}{\overset{ˆ}{U}}$consisting of only the estimated noise eigenvectors corresponding to the eigenvectors with the lowest magnitude eigenvalues.

The algorithm in this example works by iteratively minimizing the cost function relative to$A\left(\theta \right),C,\text{or}G.$

Algorithm Overview

This example implements the iterative method proposed by Friedlander and Weiss. The algorithm uses the MUSIC spatial spectrum as a starting point and improves the DOA estimates by iteratively updating the mutual coupling and calibration parameters.

Algorithm Steps:

1. Initialize - Start with nominal or previously measured values of $C$ and $G$. In this case, we initialize both to the identity matrix.

2. Estimate DOAs: Apply the MUSIC algorithm to estimate $\theta$ using the current estimates of $C$ and $G$.

3. Estimate Calibration Matrix: Estimate $G$ by minimizing $J$ keeping $a\left({\theta }_n\right)$ and $C$ constant.

4. Estimate MCM: Estimate $C$ by minimizing $J$ keeping $a\left({\theta }_n\right)$ and $G$ constant.

5. Repeat: Repeat steps 2-4 until some stopping criteria is hit.

Uniform Linear Array (ULA) Example

In this part of the example, we demonstrate the above algorithm using a ULA. We place some restrictions on the structure of the MCM in this example based on the array geometry. In reality, if the true MCM does not have such an exploitable structure, another algorithm will need to be used to form an estimate.

In the case of the ULA, we can assume that the MCM has a banded Toeplitz matrix. This is based on three underlying assumptions:

1. The assumption of symmetry, where the coupling between element i and j is the same as the coupling between the element j and i.

2. The assumption that the mutual coupling between two elements is based only on the distance between those two elements.

3. The assumption that at some distance between two elements, the mutual coupling is negligable.

If these assumptions are held to be true, it substantially reduces the number of unknown values in the MCM that we need to estimate. For example, if we have an eight-element array and assume that there is only substantial coupling between an element and two of its neighbors, we can investigate the coupling between element 4 and its neighbors:

By making the three assumptions listed above, the coupling coefficients actually just get repeated in each column of the coupling matrix. Therefore, despite $C$ being an 8x8 matrix, it has a banded Toeplitz structure that can actually be reduced to only two unknown values:

$$C=\left[\begin{array}{cccccccc}1& c_{45}& c_{46}& 0& 0& 0& 0& 0\\ c_{45}& 1& c_{45}& c_{46}& 0& 0& 0& 0\\ c_{46}& c_{45}& 1& c_{45}& c_{46}& 0& 0& 0\\ 0& c_{46}& c_{45}& 1& c_{45}& c_{46}& 0& 0\\ 0& 0& c_{46}& c_{45}& 1& c_{45}& c_{46}& 0\\ 0& 0& 0& c_{46}& c_{45}& 1& c_{45}& c_{46}\\ 0& 0& 0& 0& c_{46}& c_{45}& 1& c_{45}\\ 0& 0& 0& 0& 0& c_{46}& c_{45}& 1\end{array}\right]$$

Enforcing this structure for $C$ makes it easier to obtain a good estimate of the MCM elements, meaning that this algorithm works well for ULAs.

ULA Implementation

In this section, we implement the algorithm and investigate whether the DOA, MCM, and calibration matrix estimates are improved. In this case, we step for a certain number of iterations, although some reasonable stopping criteria can be developed such as stopping when the cost function is no longer being reduced by some threshold value.

Start by setting up the array and signal geometry. The array is an 8 element ULA with half wavelength spacing. Each receiver channel has a noise figure of 10 dB.

% Set random seed
rng(1);

% Create Array
fc = 10e9;
lambda = freq2wavelen(fc);
m = 8;
array = phased.ULA(m,lambda/2);
elPos = getElementPosition(array)/lambda;

% Create receiver
rx = phased.Receiver(NoiseFigure=10);

There are three signals arriving to the array at different azimuth angles. Note that for arrays with more antenna elements, more signals may be required to for adequate estimates.

% Create 3 signals arriving at different angles
nSnapshots = 1e5;
sigAngs = [-45 15 60];
nSig = length(sigAngs);

Generate the true mutual coupling values and calibration errors using a helper function. We assume that mutual coupling is non-zero for only the two immediately adjacent elements.

% Generate calibration errors
cal = helperGenerateCalErrors(m);

% Generate mutual coupling matrix
nNzDiag = 2;
mcm = helperGenerateUlaMcm(m,nNzDiag);

Simulate the signals arriving at the array and calculate the noise eigenvectors.

% Collect plane wave
sig = sensorsig(elPos,nSnapshots,sigAngs,db2pow(-30)).';
sig = (mcm * cal * sig).';

% Apply noise
sig = rx(sig).';

% Compute data covariance
Rx = sig*sig'/size(sig,2);
    
% Compute eigenvectors, remove signal eigenvectors leaving only noise
% eigenvectors.
[V,D] = eig(Rx);
h = diag(D);
[~,hidx] = sort(h,'descend');
U = V(:,hidx(nSig+1:end));

Initialize estimates of the calibration errors and mutual coupling matrix. In this case we use the identity matrix for both initial estimates. In reality, you may have some previously measured values.

% Initialize calibration and MCM estimates
estCal = eye(m);
estMcm = eye(m);

Initialize variables used in our loop iterations.

% Perform 30 iterations
nIter = 30;

% Check DOA from -90 to 90 in 0.01 degree increments
testAngs = -90:0.01:90;

% Initialize storage variables
doaSpectrums = cell(1,nIter);
estDoas = cell(1,nIter);
estCals = cell(1,nIter);
estMcms = cell(1,nIter);

Iteratively estimate DOA, MCM, and calibration errors.

for i = 1:nIter
    % Calculate the spatial spectrum
    sv = steervec(elPos,testAngs);
    resp = U'*estMcm*estCal*sv;
    spectrum = abs(1./sum(abs(resp).^2,1));

    % Calculate the DOA by finding the highest peaks
    [~,peaklocs] = findpeaks(spectrum,testAngs,"MinPeakDistance",5,"SortStr","descend");
    estDoa = peaklocs(1:nSig);

    % Estimate cal values holding mcm and DOA fixed
    estCal = helperGetCalEst(U,estMcm,estDoa,elPos);
    
    % Estimate mutual coupling matrix holding cal and DOA fixed
    estMcm = helperGetMcmEstUla(U,estCal,estDoa,elPos,nNzDiag);

    % Save calculated values
    doaSpectrums{i} = spectrum;
    estDoas{i} = estDoa;
    estCals{i} = estCal;
    estMcms{i} = estMcm;
end

Once we have run the algorithm, plot the results of our estimates to investigate whether the algorithm performed adequately.

helperPlotResults(sigAngs,cal,mcm,doaSpectrums,estDoas,estCals,estMcms,testAngs);

We can see that as the algorithm advances, the spectral peaks become sharper, and our estimates of the DOA, calibration matrix and MCM improve.

Uniform Circular Array (UCA) Example

We can perform the same processing steps in the case of the UCA. In this case however, the structure of the MCM matrix is slightly different than the ULA.

For the UCA, we make the same foundational assumptions as in the ULA case: the mutual coupling matrix is symmetric, the coupling strength depends only on the distance between elements, and at some distance the mutual coupling between elements is negligible. However, due to the circular geometry of the UCA, these assumptions result in a circulant mutual coupling matrix rather than a banded Toeplitz matrix. This is because a circular array has no boundaries—every element has neighbors on both sides—so each row (and column) of the mutual coupling matrix has the same pattern of non-zero entries.

For example, in an eight-element circular array where substantial coupling occurs only between each element and its two nearest neighbors, the mutual coupling pattern can be visualized as follows:

By making these assumptions, the coupling coefficients get repeated in each column of the coupling matrix. Therefore, despite $C$ being an 8x8 matrix, it has a circulant structure that can actually be reduced to only two unknown values:

$$C=\left[\begin{array}{ccccccc}1& c_{12}& c_{13}& 0& 0& c_{13}& c_{12}\\ c_{12}& 1& c_{12}& c_{13}& 0& 0& c_{13}\\ c_{13}& c_{12}& 1& c_{12}& c_{13}& 0& 0\\ 0& c_{13}& c_{12}& 1& c_{12}& c_{13}& 0\\ 0& 0& c_{13}& c_{12}& 1& c_{12}& c_{13}\\ c_{13}& 0& 0& c_{13}& c_{12}& 1& c_{12}\\ c_{12}& c_{13}& 0& 0& c_{13}& c_{12}& 1\end{array}\right]$$

Assuming this structure for $C$ makes it easier to obtain a good estimate of the MCM elements, meaning that this algorithm works well for UCAs.

UCA Implementation

In this section, we implement the algorithm for the UCA. The approach is very similar to that of the ULA with a slight variation when minimizing the cost function with respect to the MCM due to the circulant structure.

% Set random seed
rng(1);

% Create Array
fc = 10e9;
lambda = freq2wavelen(fc);
m = 8;
spacing = lambda/2;
r = (spacing/2)/sin(pi/m);
array = phased.UCA(m,r,ArrayNormal='z');
elPos = getElementPosition(array)/lambda;

% Create receiver
rx = phased.Receiver(NoiseFigure=10);

Once again there are three signals arriving to the array at different azimuth angles. Note that for arrays with more antenna elements, more signals may be required to for adequate estimates.

% Create 3 signals arriving at different angles
nSnapshots = 1e5;
sigAngs = [-45 15 60];
nSig = length(sigAngs);

Generate the true mutual coupling values and calibration errors using a helper function. We assume that mutual coupling is non-zero for only the two immediately adjacent elements.

% Generate calibration errors
cal = helperGenerateCalErrors(m);

% Generate mutual coupling matrix
nNzDiag = 2;
mcm = helperGenerateUcaMcm(m,nNzDiag);

Simulate the signals arriving at the array and calculate the noise eigenvectors.

% Collect plane wave
sig = sensorsig(elPos,nSnapshots,sigAngs,db2pow(-30)).';
sig = (mcm * cal * sig).';

% Apply noise
sig = rx(sig).';

% Compute data covariance
Rx = sig*sig'/size(sig,2);
    
% Compute eigenvectors, remove signal eigenvectors leaving only noise
% eigenvectors.
[V,D] = eig(Rx);
h = diag(D);
[~,hidx] = sort(h,'descend');
U = V(:,hidx(nSig+1:end));

Initialize estimates of the calibration errors and mutual coupling matrix. In this case we use the identity matrix for both initial estimates. In reality, you may have some previously measured values.

% Initialize calibration and MCM estimates
estCal = eye(m);
estMcm = eye(m);

Initialize variables used in our loop iterations.

% Perform 30 iterations
nIter = 30;

% Check DOA from -90 to 90 in 0.01 degree increments
testAngs = -90:0.01:90;

% Initialize storage variables
doaSpectrums = cell(1,nIter);
estDoas = cell(1,nIter);
estCals = cell(1,nIter);
estMcms = cell(1,nIter);

Iteratively estimate DOA, MCM, and calibration errors.

for i = 1:nIter
    % Calculate the spatial spectrum
    sv = steervec(elPos,testAngs);
    resp = U'*estMcm*estCal*sv;
    spectrum = abs(1./sum(abs(resp).^2,1));

    % Calculate the DOA by finding the highest peaks
    [~,peaklocs] = findpeaks(spectrum,testAngs,"MinPeakDistance",5,"SortStr","descend");
    estDoa = peaklocs(1:nSig);

    % Estimate cal values holding mcm and DOA fixed
    estCal = helperGetCalEst(U,estMcm,estDoa,elPos);
    
    % Estimate mutual coupling matrix holding cal and DOA fixed
    estMcm = helperGetMcmEstUCA(U,estCal,estDoa,elPos,nNzDiag);

    % Save calculated values
    doaSpectrums{i} = spectrum;
    estDoas{i} = estDoa;
    estCals{i} = estCal;
    estMcms{i} = estMcm;
end

Once we have run the algorithm, plot the results of our estimates to investigate whether the algorithm performed adequately.

helperPlotResults(sigAngs,cal,mcm,doaSpectrums,estDoas,estCals,estMcms,testAngs);

We can see that as the algorithm advances, the spectral peaks become sharper, and our estimates of the DOA, calibration matrix and MCM become better.

Conclusion

This example demonstrated how to perform direction finding in the presence of mutual coupling and calibration uncertainties using a self-calibrating iterative approach.

Key Takeaways:

• Mutual coupling and calibration errors can significantly degrade DOA estimation.

• When the structure of the mutual coupling matrix is known (banded Toeplitz for ULAs, circulant for UCAs), it can be exploited to estimate and correct for these effects.

• The method iteratively refines the system model using only received data, requiring no prior knowledge of source directions.

Limitations:

• Performance depends on accurate assumptions about matrix structure.

• Requires adequate signal-to-noise ratio and sufficient data snapshots.

• Convergence is not guaranteed in all scenarios and may require tuning.

This approach is especially valuable for systems where environmental or hardware conditions vary over time, necessitating online self-calibration.

References

[1] B. Friedlander and A. J. Weiss. “Direction Finding in the Presence of Mutual Coupling,” IEEE Transactions on Antennas and Propagation, vol. 39, no. 3, pp. 273–284, Mar. 1991.

[2] R. Schmidt. "Multiple emitter location and signal parameter estimation." IEEE transactions on antennas and propagation 34, no. 3 (1986): 276-280.

Helper Functions

function z = helperGenerateCalErrors(m)
    % Generate initial amplitude errors.
    ampErr = (rand(m,1)-0.5)*0.2*sqrt(12)+1;
    
    % Generate phase errors
    phaseErr = deg2rad((rand(m,1)-0.5)*20);
    
    % Generate calibration errors
    z = diag(ampErr .* exp(1i * phaseErr));
end

function mcm = helperGenerateUlaMcm(m,nNzDiag)
    % Initialize coupling coefficients (1 on diagonal, decreasing off-diagonal)
    mainDiag = 1;

    % Exponentially decreasing amplitudes
    offDiagAmplitudes = (0.4 .^ (1:nNzDiag)).';

    % Generate off diagonal phases.
    offDiagPhases = deg2rad((rand(nNzDiag,1)-0.5)*20);
    
    % Create off diagonal coefficients
    offDiagCoeffs = offDiagAmplitudes.*exp(1i*offDiagPhases);

    % Initialize a single set of coefficients, these go in one column
    coeff = [mainDiag;offDiagCoeffs];

    % Create banded toeplitz
    mcm = helperBandedToeplitz(coeff,m);
end

function toep = helperBandedToeplitz(coeff,m)
    % Assume that the first value in coeff is the center
    fullcoeff = [flip(coeff(2:end));coeff];
    nNzDiag = length(coeff)-1;
    
    % Create complex coefficients for diagonals
    toep = zeros(m);
    for k = 1:m
        coeffidxs = (k-nNzDiag):(k+nNzDiag);
        coeffkeep = coeffidxs >= 1 & coeffidxs <= m;
        toep(coeffidxs(coeffkeep),k) = fullcoeff(coeffkeep);
    end
end

function mcm = helperGenerateUcaMcm(m,nNzDiag)
    % Initialize coupling coefficients (1 on diagonal, decreasing off-diagonal)
    mainDiag = 1;

    % Exponentially decreasing amplitudes
    offDiagAmplitudes = (0.4 .^ (1:nNzDiag)).';

    % Generate off diagonal phases.
    offDiagPhases = deg2rad((rand(nNzDiag,1)-0.5)*30);
    
    % Create off diagonal coefficients
    offDiagCoeffs = offDiagAmplitudes.*exp(1i*offDiagPhases);

    % Initialize a single set of coefficients, these go in one column
    coeff = [mainDiag;offDiagCoeffs];

    % Create circulant
    mcm = helperCirculant(coeff,m);
end

function circ = helperCirculant(coeff,m)
    % Assume that the first value in coeff is the center
    nZero = m-(length(coeff)*2-1);
    col = [coeff;zeros(nZero,1);flip(coeff(2:end))];
    
    % Create circ shift matrix
    circ = zeros(m);
    for i = 0:m-1
        circ(:,i+1) = circshift(col,i);
    end
end

function calEst = helperGetCalEst(U,estMcm,estDoa,elPos)
    % Calculate calibration estimate. See [1] for additional details.
    nSig = length(estDoa);
    Z = 0;
    for n = 1:nSig
        Q = diag(steervec(elPos,estDoa(n)));
        Z = Z + Q'*estMcm'*(U*U')*estMcm*Q;
    end

    % Calculate delta
    Zinv = Z^-1;
    m = size(estMcm,1);
    w = zeros(m,1);
    w(1) = 1;
    delta = (Zinv*w)/(w'*Zinv*w);

    % Diagonalize delta to get the cal estimate
    calEst = diag(delta);
end

function estMcm = helperGetMcmEstUla(U,estCal,estDoa,elPos,nNzDiag)
    % Calculate MCM estimate. See [1] for additional details.
    nSig = length(estDoa);
    m = size(estCal,1);
    G = 0;
    for n = 1:nSig
        a = estCal*steervec(elPos,estDoa(n));
        Q = helperBuildQULA(a,m,nNzDiag);
        G = G + Q'*(U*U')*Q;
    end

    % Calculate c
    W = zeros(nNzDiag+1,1);
    W(1,1) = 1;
    u = 1;
    u(1) = 1;
    c = (G^-1)*W*((W.'*(G^-1)*W)^-1)*u;

    % Create banded toeplitz
    estMcm = helperBandedToeplitz(c,m);
end

function Q = helperBuildQULA(a,m,nNzDiag)
    % Create Q for a banded toeplitz matrix. See [1] for additional
    % details.
    Q = zeros(m,nNzDiag+1);
    for p = 1:m
        for q = 1:nNzDiag+1
            if p + q <= m + 1
                Q(p,q) = Q(p,q) + a(p+q-1);
            end

            if p >= q && q >= 2
                Q(p,q) = Q(p,q) + a(p-q+1);
            end
        end
    end
end

function estMcm = helperGetMcmEstUCA(U,estCal,estDoa,elPos,nNzDiag)
    % Calculate MCM estimate. See [1] for additional details.
    nSig = length(estDoa);
    m = size(estCal,1);
    G = 0;
    for n = 1:nSig
        a = estCal*steervec(elPos,estDoa(n));
        Q = helperBuildQUCA(a,m,nNzDiag);
        G = G + Q'*(U*U')*Q;
    end

    % Calculate c
    W = zeros(nNzDiag+1,1);
    W(1,1) = 1;
    u = 1;
    u(1) = 1;
    c = (G^-1)*W*((W.'*(G^-1)*W)^-1)*u;

    % Create banded toeplitz
    estMcm = helperCirculant(c,m);
end

function Q = helperBuildQUCA(a,m,nNzDiag)
    % Create Q for a banded toeplitz matrix. See [1] for additional
    % details.
    Q = zeros(m,nNzDiag+1);
    if mod(m/2,0)
        l = m/2;
    else
        l = (m+1)/2;
    end
    for p = 1:m
        for q = 1:nNzDiag+1
            if p + q <= m + 1
                Q(p,q) = Q(p,q) + a(p+q-1);
            end

            if p >= q && q >= 2
                Q(p,q) = Q(p,q) + a(p-q+1);
            end

            if p < q && q <= l
                Q(p,q) = Q(p,q) + a(m+1+p-q);
            end

            if 2 <= q && q <= l && (p+q) >= (m+2)
                Q(p,q) = Q(p,q) + a(p+q-m-1);
            end
        end
    end
end

function helperPlotResults(trueDoa,trueCal,trueMcm,doaSpectrums,estDoas,estCals,estMcms,testAngs)
    % Plot Spectrum
    specAx = axes(figure);
    helperPlotDoaSpectrum(testAngs,doaSpectrums,trueDoa,specAx);

    % Create figure for plotting errors
    tl = tiledlayout(figure,"vertical");

    % Plot DOA Error
    doaAx = nexttile(tl);
    helperPlotDoaError(trueDoa,estDoas,doaAx);

    % Plot Calibration Error
    calAx = nexttile(tl);
    helperPlotCalibrationEstimate(trueCal,estCals,calAx);

    % Plot MCM Error
    mcmAx = nexttile(tl);
    helperPlotMcmEstimate(trueMcm,estMcms,mcmAx);
end

function helperPlotDoaSpectrum(testAngs,doaSpectrums,trueDoa,ax)
    nSpec = length(doaSpectrums);
    toPlot = [1 5 10 20 30 50 100];
    toPlot = toPlot(toPlot<=nSpec);
    hold(ax,"on");
    for p = toPlot
        plot(ax,testAngs,mag2db(doaSpectrums{p}),DisplayName=[num2str(p)]);
    end
    xline(ax,trueDoa,LineStyle="--",Color='k',DisplayName='DoA');
    title(ax,"Spatial Spectrum After N Iterations")
    xlabel(ax,"Iteration");
    ylabel(ax,"Spectrum Magnitude (dB)");
    legend(ax,Location="eastoutside")
end

function helperPlotDoaError(trueDoa,estDoas,ax)
    % Plot the DOA error over iterations.
    nEst = length(estDoas);
    estimateVal = zeros(1,nEst);
    trueDoa = sort(trueDoa);
    for i = 1:nEst
        estDoa = sort(estDoas{i});
        estimateVal(i) = norm(trueDoa-estDoa);
    end

    % Plot goodness at each iteration
    plot(ax,estimateVal);
    xlabel(ax,"Iteration");
    ylabel(ax,"DOA Error (Degree)");
    title(ax,"DOA Estimate Error vs. Iteration");
end

function helperPlotCalibrationEstimate(trueCal,allEstCals,ax)
    % The calibration should estimate the true cal, with the first value
    % adjusted to be 1.
    trueCal = trueCal/trueCal(1,1);

    % Characterize the difference between each estimate and true by taking the
    % norm of the delta
    nEst = length(allEstCals);
    estimateVal = zeros(1,nEst);
    for i = 1:nEst
        estimateVal(i) = norm(trueCal-allEstCals{i})/norm(trueCal);
    end

    % Plot goodness at each iteration
    plot(ax,estimateVal);
    xlabel(ax,"Iteration");
    ylabel(ax,"Calibration Error");
    title(ax,"Calibration Estimate Error vs. Iteration");
end
   
function helperPlotMcmEstimate(trueMcm,allEstMcm,ax)
    % Characterize the difference between each estimate and true by taking the
    % norm of the delta
    nEst = length(allEstMcm);
    estimateVal = zeros(1,nEst);
    for i = 1:nEst
        estimateVal(i) = norm(trueMcm-allEstMcm{i})/norm(trueMcm);
    end

    % Plot goodness at each iteration
    plot(ax,estimateVal);
    xlabel(ax,"Iteration");
    ylabel(ax,"MCM Error");
    title(ax,"MCM Estimate Error vs. Iteration");
end