clear; clc; close all;
LINEWIDTH = 3;
XLIM = [-20, 20];
%% ========== Simulation Parameters ==========
% ==================== SNR ====================
SNRdB = 0;
% ==================== Snapshots ====================
SNAPSHOTS = 1000;
% ==================== Number of sources ====================
DD = 4;
% ==================== Source distribution ====================
theta_bar = (0:DD-1).'*0.1;
%theta_bar = zeros(DD,1);
%theta_bar(:,1)=[10 20 30 40];
amp = ones(DD, 1);
% ==================== Noise power ====================
noise_std = 10^(-SNRdB/20) / sqrt(DD);
% ==================== Mutual coupling parameters ====================
c1 = 0.5*exp(0.25i*pi);
c2 = 0.5*exp(0.7i*pi)/2;
c3 = 0.5*exp(0.7i*pi)/3;
B = 3;
%% ========== Compute difference coarrays and weights ==========
TYPES = 4;
for type = 1 : TYPES
switch (type)
case 1
% ULA
S = [0; 1; 2; 3; 4; 5];
name = 'ULA';
case 2
% MRA
S = [0; 1; 6; 9; 11; 13];
name = 'MRA';
case 3
% Nested
S = [1; 2; 3; 4; 8; 12];
name = 'Nested';
case 4
% Coprime
S = [0; 2; 4; 3; 6; 9];
name = 'Coprime';
end
[V, w_on_V] = weight_function( S, 'V' );
[U, w_on_U] = weight_function( S, 'U' );
LEN_S = length(S);
LEN_U = length(U);
% ========== Weight function ==========
figure;
stem(V, w_on_V, 'LineWidth', LINEWIDTH)
xlabel('Coarray location $m$', 'interpreter', 'latex')
ylabel('$w(m)$', 'interpreter', 'latex')
xlim(XLIM); grid on; set(gca, 'YMinorGrid', 'off');
title(name, 'interpreter', 'latex')
% ========== MUSIC spectra ==========
[n1, n2] = ndgrid(S);
n1_n2_mat = n1 - n2;
n1_n2_vec = n1_n2_mat(:);
% Steering vector
V = exp(2i * pi * S * theta_bar');
% Source Process
Source = diag(amp) * (randn(DD, SNAPSHOTS) + 1i * randn(DD, SNAPSHOTS)) / sqrt(2);
% Noise Process
Noise = diag(noise_std) * (randn(LEN_S, SNAPSHOTS) + 1i * randn(LEN_S, SNAPSHOTS)) / sqrt(2);
% Coupling matrix
C_Q = eye(LEN_S);
C_Q( abs(n1_n2_mat) == 1) = c1;
C_Q( abs(n1_n2_mat) == 2) = c2;
C_Q( abs(n1_n2_mat) == 3) = c3;
% Numerical values of the matrix C
disp(['abs(C) for ', name, ':'])
disp( num2str( abs(C_Q) ) )
disp('========================================')
X_Q = C_Q * V * Source + Noise;
R_Q_SxS = X_Q * X_Q' / SNAPSHOTS;
if (type == 1)
x_D = [R_Q_SxS(1, end:-1:1).'; R_Q_SxS(2:end, 1)];
% Plot MUSIC spectrum
[P_Q, theta_bar_grid] = spectral_MUSIC(x_D, DD);
figure;
semilogy( theta_bar_grid, P_Q, 'LineWidth', LINEWIDTH)
grid on;
set(gca, 'XTick', theta_bar, 'YMinorGrid', 'off');
xlabel('$\bar{\theta}$', 'interpreter', 'latex')
ylabel('$P(\bar{\theta})$', 'interpreter', 'latex')
theta_bar_est = root_MUSIC( U, x_D, DD);
title(['ULA, Error = ', num2str(norm( theta_bar_est - theta_bar ) / sqrt(DD) )], 'interpreter', 'latex')
else
x_SxS = R_Q_SxS(:);
x_U = zeros(LEN_U, 1);
for kk = 1 : LEN_U
x_U(kk) = mean( x_SxS(n1_n2_vec == U(kk) ) );
end
% Plot MUSIC spectrum
[P_Q, theta_bar_grid] = spectral_MUSIC(x_U, DD);
figure;
semilogy( theta_bar_grid, P_Q, 'LineWidth', LINEWIDTH)
grid on;
set(gca, 'XTick', theta_bar, 'YMinorGrid', 'off');
xlabel('$\bar{\theta}$', 'interpreter', 'latex')
ylabel('$P(\bar{\theta})$', 'interpreter', 'latex')
theta_bar_est = root_MUSIC( U, x_U, DD);
title([name, ', Error = ', num2str(norm( theta_bar_est - theta_bar ) / sqrt(DD) )], 'interpreter', 'latex')
end
end
