%基于高階累計量的MUSIC算法
clear;
clc;
tic;
j = sqrt(-1);
K = 1;
M = 5;
d = 0.5;
t = [0:1/200:1]';
N = length(t);
S = [sin(2*pi*20*t)]';
theta_s = [10];
SNR = 10;
p_noise = (1/K*2)/10^(SNR/10);
n = p_noise/2*randn(M,N);
for k = 1:K
A(:,k) = exp(-j*(0:M-1)'*2*d*pi*sin(theta_s(k)*pi/180));
end
X = A*S+n;
Rcum=(kron(X,conj(X))*(kron(X,conj(X))'))/N-kron(X,conj(X))/N*kron(X,conj(X))'/N-kron(X*X'/N,conj(X*X'/N));
[V D] = svd(Rcum);
En = V(:,2:M^2); %噪聲子空間
delt_theta = 1;
theta = -90:delt_theta:90;
for num = 1:length(theta)
at = exp(-j*(0:(M-1))'*2*d*pi*sin(theta(num)*pi/180));
bt = kron(at,conj(at));
y = En'*bt;
P_index(num) = y'*y; %連續改變theta進行搜索
end
P_index = 1./P_index;
P_index = 10*log10(P_index);
plot(theta,P_index,'b');
xlabel('Angle(degree)');ylabel('Spectrum(dB)');
hold on;
grid on;
toc;