Multiplicative Array Processing

I will give a short answer first, then summarize the mathematical basis of the paper, and finally provide a MATLAB code example.

Short answer

MAP appears to be taking advantage of the same thing that root-MUSIC does, which is that in the case of a uniform linear array (ULA), the relationship between the received signals across channels has a polynomial structure. MAP uses this polynomial structure to synthesize a beampattern from a few elements as if you had more elements. This polynomial structure also explains why the nulls in an antenna pattern are much narrower than the beamwidth: The nulls are the roots of the polynomial.

I think the fundamental thing MAP is taking advantage of is that we are assuming the signal is a sine wave (i.e., narrowband), which means the phase response across a ULA is linear. Because of this, if we measure the phase difference between two adjacent array elements, that tells us the slope of the linear phase ramp across the whole array (if you aren't super familiar with phased array beamforming, this answer should help). All we need to compute the equation of a line is 2 points.

So we are using our knowledge/assumptions about the signal to improve our spectral estimate (root-MUSIC is taking advantage of the same thing in a different way). This is the heart of any super resolution technique: using knowledge of the signal we are measuring to get better resolution than the measurements alone would dictate.

However, there is a catch, which is that if you don't actually measure the signal on additional array elements, and instead assume that the observed phase difference between 2 elements holds for the rest of the array, you will be vulnerable to noise. An N-channel array allows you to improve SNR by a factor of N. You lose that advantage if you use a 2 element array and MAP instead.

The math

For a signal arriving from direction $\theta_0$ relative to array boresight, if you beamform $N$ receive channels (spaced by distance $d$) for a set of angles $\theta$, the effect on the arriving signal (ignoring amplitude) as a function of $\theta$ will be:

$$ r(\theta) = \sum_{n=0}^{N-1} e^{-j \frac{2\pi}{\lambda}nd(\sin\theta - sin\theta_0)} = \sum_{n=0}^{N-1} X^n = 1 + X + X^2 + \cdots X^{N-1} $$ where $X = e^{-j\frac{2\pi}{\lambda} d \cdot (\sin \theta - \sin \theta_0) } $ and $\lambda$ is the signal wavelength ($\frac{2\pi}{\lambda} = k_0$ in the paper).

The above demonstrates that for a ULA, computing a DoA spectrum via beamforming to directions $\theta$ amounts to evaluating a polynomial of degree $N-1$. The trick is that instead of summing powers of $X$ to evaluate the polynomial, we could multiply factors of the polynomial and get the same answer:

$$ r(\theta) = 1 + X + X^2 + \cdots X^{N-1} = (X + x_1)(X + x_2) \cdots (X + x_{N-1}) $$ where $x_n$ are the $N-1$ complex roots of the polynomial. You could also combine some of the factors and have fewer terms.

Example

Suppose $N = 4$. Then $r(\theta) = 1 + X + X^2 + X^3 = (1 + X)(1 + X^2)$

What each of these factors represents is a linear combination of the signals received on a pair of elements. Thus, the first factor $(1 + X)$ is the combination of the first element (the phase reference) and the second element. The second factor is the combination of the first and third elements.

Therefore, the factored version of the polynomial tells us that we can get the same DoA pattern by:

1. Summing the response on the first and second elements
2. Summing the response on the first and third elements
3. Multiplying the two sums

The DoA pattern $r(\theta)$ will be the same, but notice that we have only used 3 of the 4 array elements to compute it! This is the basic principle of MAP. The figure below illustrates this for a signal DoA of 30 degrees.

We can extend this to use only two elements to compute the pattern from a 4 element array by factoring the second polynomial into complex roots: $r(\theta) = (1 + X)(1 + X^2) = (1 + X)(1 + jX)(1 - jX)$. Each of these factors represents different linear combinations of the response on the first 2 elements of the array, which when multiplied together give the same pattern as the linear 4 element array.

MATLAB Code

%Scenario parameters 
c = 343;        %Propagation speed (m/s)
N = 4;          %Number of elements in array
fc = 2e3;       %Carrier frequency (Hz)
lambda = c/fc;  %Wavelength
d = lambda/2;   %Element spacing
theta_0 = 30;   %Signal direction of arrival (relative to array boresight)

%Signal channel phase 
s = exp(1j*2*pi/lambda*d*[0:N-1].'*sind(theta_0));

%Steering vectors for directions theta
thetas = -90:0.1:90;
steerVecs = exp(1j*2*pi/lambda*d*[0:N-1].'.*sind(thetas));  

%Apply channel weights (steering vector coefficients) to received signal
s_weight = conj(s).*steerVecs; 

%Compute the DoA spectrum for all theta using all N channels
rThetaN = s'*steerVecs;     %rThetaN = sum(s_weight, 1)
figure; 
plot(thetas, 20*log10(abs(rThetaN)), 'linewidth', 2); 
ylim([20*log10(N) - 50, 20*log10(N)])
xlabel('Angle of Arrival')
ylabel('Gain (dB)')
xticks(-90:30:90)

%Compute the DoA spectrum for all theta using 3 channel MAP
rTheta3 = (s_weight(1,:) + s_weight(2,:)).*(s_weight(1,:) + s_weight(3,:));
hold on; plot(thetas, 20*log10(abs(rTheta3)), 'r--', 'linewidth', 2)
legend('Linear (N = 4)', 'MAP (N = 3, PG = 4/3)', 'location', 'best')

%Compute the DoA spectrum using 2 channel MAP
rTheta2 = (s_weight(1,:) + s_weight(2,:)).*(s_weight(1,:) + 1j*s_weight(2,:)).*(s_weight(1,:) - 1j*s_weight(2,:));

figure; 
plot(thetas, 20*log10(abs(rThetaN)), 'linewidth', 2); 
ylim([20*log10(N) - 50, 20*log10(N)])
xlabel('Angle of Arrival')
ylabel('Gain (dB)')
hold on; plot(thetas, 20*log10(abs(rTheta2)), 'g--', 'linewidth', 2)
legend('Linear (N = 4)', 'MAP (N = 2, PG = 4/2)', 'location', 'best')
xticks(-90:30:90)