1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
%% MUSIC算法---空间谱估计
%----单信源空间谱估计
clear all;
close all;
clc;

%% 参数设置
c=3e8;              % 光速
fc=500e6;           % 载波频率
fs=1e9;             % 采样频率1GHz
SNR=10;             % 信噪比dB
M=16;               % 阵元个数
lamda=c/fc;         % 波长
d=lamda/2;          % 阵元间距
COUNT=1024;         % 快拍数
n=[0:M-1]';
thetas=[-30 35 60];     % 波束指向
f_dp=[300 400 500]*1e6; % 与波束指向对应的信号频率
Num_source=length(thetas);          % 信源个数

bujing=0.1;
theta=-90:bujing:90;                % 空间扫描角度范围
t=(0:COUNT-1)/fs;

P_noise=1;              % 噪声功率
P_signal=10^(SNR/10);   % 信号功率

%% 得到输入信号的协方差矩阵

for p=1:length(thetas)
    Rt(p,:)=P_signal*exp(1i*2*pi*f_dp(p)*t);
end
Vs=exp(1j*2*pi*fc*n*d*sind(thetas)/c);                   % 方向矢量
noise=(rand(M,COUNT)+1j*rand(M,COUNT))/sqrt(2);
X=Vs*Rt+noise;
R_xx=X*X'/COUNT;         % 阵列接收信号的协方差矩阵


%% MUSIC 算法
[V,D]=eig(R_xx);         % 特征值D和特征向量V;
U=fliplr(V);             % 将矩阵A的列绕垂直轴进行左右翻转
for p=1:length(theta)
    a=exp(1j*2*pi/lamda*n*d*sind(theta(p))); % 方向矢量
    G_n=U(:,(Num_source+1:M) );       % 噪声子空间的特征向量构成的矩阵
    P_temp(p)=1./(a'*G_n*G_n'*a);     % MUSIC算法空间谱
end

P=abs(P_temp);
PP=10*log10(P/max(P)); % 归一化功率
figure(1)
plot(theta,PP);grid on;
xlabel('入射角度/度');ylabel('信号功率/dB');
title('MUSIC谱');