非圆信号DOA估计(NC-MUSIC附上源码,emmm,发现大家都喜欢有code)

EstelleYuan

已于 2023-04-07 20:39:04 修改

阅读量1.2k

点赞数 23

分类专栏： DOA 非圆信号文章标签：算法

于 2021-12-10 22:48:48 首次发布

本文链接：https://blog.csdn.net/jianai223/article/details/121866335

版权

DOA 同时被 2 个专栏收录

6 篇文章 6 订阅

订阅专栏

非圆信号

2 篇文章 1 订阅

订阅专栏

本文探讨了NC-MUSIC算法在信噪比较低环境下对多目标方向估计的准确性问题。通过MATLAB代码展示了算法执行过程，并指出其在SNR=0dB时估测偏差大，与MUSIC算法相比抗噪性能较差。作者寻求解释这一现象并邀请读者分享见解。

摘要由CSDN通过智能技术生成

% code by Estelle Yuan
clear all;
close all;
clc;
derad = pi/180;
radeg = 180/pi;
twpi = 2 * pi;
M = 3;  % ULA阵元数

lamda = 1;
dd = lamda / 2;  %阵元间距
d = 0 : dd : (M - 1) * dd;

theta = [10 30 50 70];
phi = [20 40 60 80];
% theta = [-0.5 0.5];
K = length(theta);  % 信元数
A = exp(1j * d.' * twpi * sin(theta * pi / 180) );    % 阵列流型

%%
Nsub = 300;   % 快拍数
s1 = randn(1, Nsub);
s1 = s1 - mean(s1).* ones(1, Nsub);
s2 = randn(1, Nsub);
s2 = s2 - mean(s2).* ones(1, Nsub);
s3 = randn(1, Nsub);
s3 = s3 - mean(s3).* ones(1, Nsub);
s4 = randn(1, Nsub);
s4 = s4 - mean(s4).* ones(1, Nsub);
ns1 = s1 .* exp(1j * deg2rad(phi(1)));
ns2 = s2 .* exp(1j * deg2rad(phi(2)));
ns3 = s3 .* exp(1j * deg2rad(phi(3)));
ns4 = s4 .* exp(1j * deg2rad(phi(4)));
ns = [ns1;ns2;ns3;ns4];
x = A * ns;
snr = 5;
x = awgn(x, snr);
y = [x;conj(x)];
Ry = y * y' / Nsub;

[EVx, Dx] = eig(Ry);                                  % 特征值分解 MUSIC-LIKE 
EVAx = diag(Dx)';                                       
[EVAx, I1x] = sort(EVAx);                                  % 特征值从小到大排列
EVAx = fliplr(EVAx);                                     % 左右翻转，从大到小排列
EVx = fliplr(EVx(:, I1x));                                     % 对应特征矢量排序
En1 = EVx(1 : M, K + 1 : 2 * M);
En2 = EVx(M + 1 : 2 * M, K + 1 : 2 * M);

abs(conj(En2)*En2.'-En1*En1')
searching_doa = -90 : 0.1 : 90;%线阵的搜索范围为-90~90度
for ii = 1:length(searching_doa)
    a_theta = exp(1j * (0 : M - 1)' * 2 * pi * dd * sin(pi * searching_doa(ii) / 180));
    f1 = a_theta' * En1 * En1' * a_theta;
    f2 = a_theta.' * En2 * En2' * a_theta;
    f3 = a_theta' * En1 * En2' * conj(a_theta);
    f4 = a_theta' * En2 * En2' * a_theta;
    Pmusic1(ii, 1) = f1 + f2 - 2 * abs(f3);
    Pmusic2(ii, 1) = f1 - abs(f3);
    Pmusic3(ii, 1) = f1 * f1' - f2 * f4;
end

% figure(2);
% plot(searching_doa, abs(Pmusic2),'r','linewidth',2);
figure(21);
P_MUSIC_dB2 = 10 * log(abs(Pmusic2));
plot(searching_doa, P_MUSIC_dB2,'r','linewidth',2);

[val, ind] = findpeaks(-P_MUSIC_dB2(:, 1));
[~, id] = sort(val, 'descend');
idd = ind(id);
doa = sort(searching_doa(idd(1 : K)))

doa =

10.1000 32.0000 49.7000 67.2000