我正在努力在 MATLAB 中模拟 MUSIC 算法以确定模拟信号的方位角到达。我的虚拟天线阵列由排列在一个平面上的 4 个天线组成。

虽然我的算法“工作”，因为 MUSIC 在很大程度上选择了对到达角的良好估计，但我注意到我的相关矩阵比算法声称的等级更高。这导致少一个非零特征值，因此对应于噪声子空间的特征向量少一个。此外，这个额外的非零特征值并不接近于零，并且将相应的特征向量作为噪声子空间的一部分会大大降低 MUSIC 的性能。

当我查看 2D 案例时，我注意到了同样的行为；即，在搜索方位角和仰角时。我在下面复制了我的代码。谁能阐明这里发生了什么？

% Assume 0 degrees due north, with positive clockwise orientation
% Positive phase difference indicates that antennna 1 recieves the waveform
% before the respective antenna (2,3, or 4)
% 
% Fix antenna 1 at the origin. From the x-axis, rotate antenna 2 60 degrees,
% rotate antenna 3 90 degrees, and rotate antenna 4 120 degrees. The
% optimal antenna spacing is still to be determined.
% 
%                         | 0 deg         --> +
%                         |
%                         |
%                         |
%             ---------- A1 ----------
%                   -60 ' | `  60
%                      '  |  `
%                     '   A3  `
%                  A4'         `A2


num_antennas = 4;
num_sig_sources = 1;
RF_freq = 600e6;        % freq of incoming signal
lambda = 3e8/RF_freq;   % wavelength of incoming signal
fs = 100e6;             % sampling frequency
samples = 200;          % number of snapshots
cos_freq = 1e6;         % freq to downconvert to

ant_spacings = [.25, .98, .91]; % antenna spacing

% generate_lookup_table just calculates a set of theoretical phase
% differences corresponding to some random AoA for the purpose of signal
% generation.
n = randi([0,360],1);
[AoAs, PDs1to2, PDs1to3, PDs1to4] = generate_lookup_table(-90, 90, .5, lambda, ant_spacings(1), ant_spacings(2), ant_spacings(3));
AoA = AoAs(n);

% generate signal
t = linspace(0,samples/fs,samples);
recived_signal_ant_1 = cos(2*pi.*t.*cos_freq);
recived_signal_ant_2 = cos(2*pi.*t.*cos_freq + PDs1to2(n));
recived_signal_ant_3 = cos(2*pi.*t.*cos_freq + PDs1to3(n));
recived_signal_ant_4 = cos(2*pi.*t.*cos_freq + PDs1to4(n));


% build sample correlation matrix
cov_mat = zeros(4,4);
for i = 1:samples
    X = [recived_signal_ant_1(i);
         recived_signal_ant_2(i); 
         recived_signal_ant_3(i);
         recived_signal_ant_4(i)];
    cov_mat = cov_mat + X*X';
end
cov_mat = cov_mat*(1/samples);

% grab eigenvalues and eigenvectors
[eig_vecs, eig_vals] = eig(cov_mat);

% generate matrix w/ noise eigenvector columns. An eigenvector corresponds
% to noise if the associated eigenvalue is close to 0. 
col_num = 1;
for i = 1:4
    if eig_vals(i,i) < 1e-4
        noise_vecs(:,col_num) = eig_vecs(:,i);
        col_num = col_num + 1;
    end
end

% walk t
theta = linspace(-90,90,200);
spectrum_func = zeros(1,100);
for i = 1:length(theta)

    % steering vector - only taking real part as of now (Re(e^i*phasediff) = cos(phase_diff))
    steering_vec = [1;
        cos((-2*pi*ant_spacings(1)/lambda)*sin(theta(i)*(pi/180) - pi/3));
        cos((-2*pi*ant_spacings(2)/lambda)*sin(theta(i)*(pi/180) - pi/2));
        cos((-2*pi*ant_spacings(3)/lambda)*sin(theta(i)*(pi/180) - 2*pi/3))];
    spectrum_func(i) = 1/(((steering_vec'*noise_vecs)*noise_vecs')*steering_vec);
end

% grab theta corresponding to highest peak, round to nearest half integer
[mx,argmax] = max(spectrum_func);
estimated_AoA = floor(2*theta(argmax)+.5)/2;
fprintf("Estimated AoA: %.1f \t Hard-coded AoA: %.1f.\n", estimated_AoA, AoA)

plot(theta,spectrum_func)
title('MUSIC Psuedospectrum')
xlabel('Angle of Arrival')

编辑：我添加了 generate_lookup_table() 的源代码：

function [AoAs, PDs1to2, PDs1to3, PDs1to4] = generate_lookup_table(start, stop, step, lambda, d1to2, d1to3, d1to4)
% d1to2 = spacing between antennas 1 and 2

num_pts = (stop-start)/step + 1;
AoAs = zeros(1,num_pts);
PDs1to2 = zeros(1,num_pts);
PDs1to3 = zeros(1,num_pts);
PDs1to4 = zeros(1,num_pts);

for i = 1:num_pts 
    AoAs(i) = start + (i-1)*step;

    % phase differences
    value_to_push_to_2 = (-2*pi*d1to2/lambda)*sin(AoAs(i)*(pi/180) - pi/3);
    value_to_push_to_3 = (-2*pi*d1to3/lambda)*sin(AoAs(i)*(pi/180) - pi/2);
    value_to_push_to_4 = (-2*pi*d1to4/lambda)*sin(AoAs(i)*(pi/180) - 2*pi/3);

    % increment/decrement until in -180 to 180
    if value_to_push_to_2 > 2*pi
        n = ceil(value_to_push_to_2/(2*pi) - 1);
        value_to_push_to_2 = value_to_push_to_2 - 2*pi*n;
    end

    if value_to_push_to_2 > pi
        value_to_push_to_2 = value_to_push_to_2 - 2*pi;
    end

    if value_to_push_to_3 > 2*pi
        n = ceil(value_to_push_to_3/(2*pi) - 1);
        value_to_push_to_3 = value_to_push_to_3 - 2*pi*n;
    end

    if value_to_push_to_3 > pi
        value_to_push_to_3 = value_to_push_to_3 - 2*pi;
    end

    if value_to_push_to_4 > 2*pi
        n = ceil(value_to_push_to_4/(2*pi) - 1);
        value_to_push_to_4 = value_to_push_to_4 - 2*pi*n;
    end

    if value_to_push_to_4 > pi
        value_to_push_to_4 = value_to_push_to_4 - 2*pi;
    end

    if value_to_push_to_2 < -2*pi
        n = ceil(-1 - value_to_push_to_2/(2*pi));
        value_to_push_to_2 = value_to_push_to_2 + 2*pi*n;
    end

    if value_to_push_to_2 < -pi
        value_to_push_to_2 = value_to_push_to_2 + 2*pi;
    end

    if value_to_push_to_3 < -2*pi
        n = ceil(-1 - value_to_push_to_3/(2*pi));
        value_to_push_to_3 = value_to_push_to_3 + 2*pi*n;
    end

    if value_to_push_to_3 < -pi
        value_to_push_to_3 = value_to_push_to_3 + 2*pi;
    end
    if value_to_push_to_4 < -2*pi
        n = ceil(-1 - value_to_push_to_4/(2*pi));
        value_to_push_to_4 = value_to_push_to_4 + 2*pi*n;
    end

    if value_to_push_to_4 < -pi
        value_to_push_to_4 = value_to_push_to_4 + 2*pi;
    end

    if value_to_push_to_4 < -2*pi
        n = ceil(-1 - value_to_push_to_4/(2*pi));
        value_to_push_to_4 = value_to_push_to_4 + 2*pi*n;
    end

    if value_to_push_to_4 < -pi
        value_to_push_to_4 = value_to_push_to_4 + 2*pi;
    end

    PDs1to2(i) = value_to_push_to_2;
    PDs1to3(i) = value_to_push_to_3;
    PDs1to4(i) = value_to_push_to_4;
end
end
```