Multidimensional harmonic retrieval¶

Data¶

Three sources impinge on a URA with azimuth angles of 10°, 30° and 70°, respectively, and with elevation angles of 20°, 30° and 40°, respectively. We observe 15 time samples, such that a tensor \(\ten{T}\in\C^{10\times 10\times 15}\) is obtained with \(t_{ijk}\) the observed signal from antenna \((i,j)\) sampled at time instance \(k\). Each source contributes a rank-1 term to the tensor. The vectors in the first and second mode are Vandermonde and the third mode contains the respective source signals multiplied by attenuation factors. Hence, the factor matrices in the first and second mode, denoted as \(\mat{A}\) and \(\mat{E}\), are Vandermonde matrices:

N = 15; R = 3; I = 10;
azimuth_angles = [10 30 70];    % azimuth angles in degrees
elevation_angles = [20 30 40];  % elevation angles in degrees
lambda = 0.25;                  % wavelength of signal carriers
dx = 0.1;                       % sensor distance
for r = 1:R
  zar = exp(2*pi/lambda*sin(azimuth_angles(r)*pi/180)*dx*1i);
  for i = 1:I, A(i,r) = zar^(i-1); end
  zer = exp(2*pi/lambda*sin(elevation_angles(r)*pi/180)*dx*1i);
  for i = 1:I, E(i,r) = zer^(i-1); end
end

The factor matrix in the third mode is the matrix containing the attenuated sources, denoted by \(\mat{S}\):

S = round(rand(N,R))*2-1;
attenuation_vector = [1 .5 .25];
for r = 1:R
  S(:,r) = attenuation_vector(r)*S(:,r);
end

The source signals are shown in Figure 65. The generated tensor is perturbed by Gaussian noise with a signal-to-noise-ratio of 0 dB:

T = cpdgen({A,E,S});
SNR = 0;
Tnoisy = noisy(T,SNR);

Figure 66 shows two observed signals with and without noise.

Figure 65: The three BPSK sources.

Figure 66: Two observed signals with and without noise, for antennas (4,5) and (5,6).

Signal separation and direction-of-arrival estimation¶

The sources are separated by means of a basic CPD, without using the Vandermonde structure.

Uest = cpd(Tnoisy,R);
Aest = Uest{1};
Eest = Uest{2};
Sest = Uest{3};

The relative errors on the estimates of the factor matrices can be calculated with cpderr. They are 0.3127, 0.2540 and 0.1698, respectively. The cpderr command also returns estimates of the permutation matrix and scaling matrices, which can be used to fix the indeterminacies:

[err,P,D,Uestfix] = cpderr({A,E,S},Uest);
Aestfix = Uestfix{1}; errA = err(1)
Eestfix = Uestfix{2}; errE = err(2)
Sestfix = Uestfix{3}; errS = err(3)

The source signals and their recovered versions are compared in Figure 67.

Figure 67: The original and recovered source signals.

The direction-of-arrival angles can be determined using the shift invariance property of the individual Vandermonde vectors. This gives relative errors for the azimuth angles of 0.0390, 0.0362 and 0.1148, and for the elevation angles of 0.0073, 0.0067 and 0.0307:

for r = 1:R
  azimuth_z_est = Aestfix(1:end-1,i)\Aestfix(2:end,i);
  azimuth_angles_est(i) = asin(angle(azimuth_z_est)*lambda/(2*pi*dx))*(180/pi);
  elevation_z_est = Eestfix(1:end-1,i)\Eestfix(2:end,i);
  elevation_angles_est(i) = asin(angle(elevation_z_est)*lambda/(2*pi*dx))*(180/pi);
end

azimuth_err = abs(sort(azimuth_angles)-sort(azimuth_angles_est))./ ...
abs(sort(azimuth_angles))
elevation_err = abs(sort(elevation_angles)-sort(elevation_angles_est))./ ...
abs(sort(elevation_angles))

Note that one can additionally impose the Vandermonde structure in the structured data fusion framework, as illustrated in other demos such as the demo Using advanced Tensorlab features for ICA.

Missing values due to broken sensors¶

If one or more antennas are broken, the tensor \(\ten{T}\) is incomplete. Tensorlab is able to process full, sparse and incomplete tensors. The missing entries can be indicated by NaN values. We consider the equivalent of a deactivated sensor, a sensor that breaks half way the experiment, and a sensor that starts to work half way the experiment:

Tnoisy(2,3,:) = NaN;
Tnoisy(5,1,1:floor(end/2)) = NaN;
Tnoisy(9,7,ceil(end/2):end) = NaN;

The incomplete tensor is visualized in Figure 68. When computing the CPD of the incomplete tensor, relative errors of 0.3180, 0.2490 and 0.1736 on the azimuth factor matrix, the elevation factor matrix and the attenuated source matrix are obtained, respectively.

Uest = cpd(Tnoisy,R);
[err,P,D] = cpderr({A,E,S},Uest);
errA = err(1)
errE = err(2)
errS = err(3)

Figure 68: Visualization of the data tensor in the case of broken sensors.