空间平滑MUSIC算法估计相干源DOA

% coherent DOA estimation
% using classic MUSIC 
% d-D searching MUSIC d-D searching using big eigen vector 1-D searching 
% collus wang 2009-11-28

clc
clear
close all

M = 8;              % 阵元数
theta = [70100]*pi/180;       % 信号入射方向
w = 2*pi*[3030]*1e3;          % 入射信号频率
Coherent_On = 1;
d = length（theta）;                  % 信号个数
T = 2e6^（-1）;                           % 采用间隔
N = 5000;                                    % 快拍数

f = 2.4e9;                          % 载波频率
c = 3e8;                            % 光速
D = c/f/2;                          % 阵元间距
SNR_dB = 10;                % SNR in dB

A = zeros（Md）;             % 导向矩阵
for index = 1 : d
    A（: index） = [exp（ 1j*2*pi*f*[0:M-1]*D*cos（theta（index））/c ）].‘;
end

s = zeros（d N）;            % 信号
for index = 1 : d
    s（index:） = exp（ 1j*w（index）*[0:N-1]*T ）;%+ rand（1N）*2*pi
end

rho = ones（1d）;            % 相干源相关系数
rho（2:d） = exp（1j*rand（1d-1）*2*pi）;
if （Coherent_On）
    for index = 2:d
        s（index:） =  s（1 :）.*rho（index）;
    end
end

sigma_noise = 1/sqrt（10^（SNR_dB/10））;       % 加噪 以单个信号的信噪比进行计算
X = zeros（M N）;
X = A*s + sigma_noise*（randn（MN）+1j*randn（MN））/sqrt（2）;       % 接收信号

% ------------------- Beginning of Algorithms -------------
R_xx = X*X‘/N;                  % calution of R_xx
[UDia]=eig（R_xx）;              % eigen value decomposition
Dia

sum = trace（abs（Dia））;          % find big eigen values
temp = 0;
pct = 0.90                                  % if the sum of biggest eigen values exceed pct of trace we say we get the biggest values.
for index  = M : -1 :  1
    temp = temp + abs（Dia（indexindex））;
    if temp > pct*sum
        break
    end
end
sig_index = index : M
noise_index = 1: index-1;

res = 0.5;                                      % resolution of peak searching （in deg）
theta_search = 0: res*pi/180: pi;

% ------------- classic music---------------
P_music_classic = zeros（1 length（theta_search） ）;
P = U（: noise_index）*U（:noise_index）‘;
for index = 1:length（theta_search）
    E = [exp（1j*2*pi*f*[0:M-1]*D*cos（theta_search（index））/c ）].‘;
    P_music_classic（index） = 1/（E‘*P*E）;
end

figure（10）
plot（theta_search/pi*180 abs（P_music_classic））
hold on
st