MUSIC和ESPRIT算法

clear all;
close all;
%% MUSIC algorithm
% 2019-11-02
%%
dread = pi/180;
N = 3;%阵元个数
M = 2;%信源个数
K = 1024;
snr = 10;
theta =[-30,10,20];%待估角度
lamda = 1;%波长
d = 0:lamda/2:(N-1)*lamda/2;%阵列
A = exp(-1i*2*pi*d.'*sin(dread *theta));
%%
S = randn(M,K);
X = A*S;
X1 = awgn(X,snr,'measured');
Rxx = X1 * X1'/K;
% Rxx = [3 0 -2;0 3 0;-2 0 3];
[EV,D] = eig(Rxx);
EVA = diag(D)';
[EVA,I] = sort(EVA);
EV = fliplr(EV(:,I));
%%
for iang = 1:361
   angle(iang) = (iang-181)/2;
    phim =dread * angle(iang);
    a = exp(-1i * 2*pi*d*sin(phim)).';% .' is 数组转置
    En = EV(:,M+1:N);
    SP(iang) = 1/(a'* En * En'* a);  
end
%%
SP = abs(SP);
SP = 10*log10(SP);
h = plot(angle,SP);
set(h,'Linewidth',0.5);
xlabel('入射角/°');
ylabel('空间谱/dB');
set(gca,'XTick',[-100:20:100]);
grid on;
%%

%% 此算法是参照博客改动的 在理解了MUSIC算法之后，自己改动一些

clear all;
close all;
% root MUSIC algorithm
% 2019-11-02
% code by wk
dread = pi/180;
N = 10;%阵元个数
M = 3;%信源个数
K = 1024;
snr = 10;
theta =[-30,10,20];%待估角度
lamda = 1;%波长
d = 0:lamda/2:(N-1)*lamda/2;%阵列
A = exp(-1i*2*pi*d.'*sin(dread *theta));
S = randn(M,K);
X = A*S;
X1 = awgn(X,snr,'measured');
Rxx = X1 * X1'/K;
% Rxx = [3 0 -2;0 3 0;-2 0 3];
[EV,D] = eig(Rxx);
EVA = diag(D)';
[EVA,I] = sort(EVA);
EV = fliplr(EV(:,I));
En = EV(:,M+1:N);
G =En *En';
m = zeros(1,2*N-1);
for i = -(N-1):(N-1)
    m(N+i) = sum(diag(G,i));  %这一步跟自己推导的不一致？ 
end
a = roots(m);
a1 = a(abs(a)<1);
[lam,I] = sort(abs(abs(a1)-1));
f = a1(I(1:M));
source_doa=zeros(1,M);
for j = 1:M
  source_doa(j)=asin(angle(f(j))/pi)*180/pi;  
end
disp('doa');
disp(source_doa);

%%根据博客看明白了公式推导 然后上面提出自己的疑问，上面在于b1n和bn1（多项式系数）和自己推导不同，按照自己的，%结果为正确的相反数。

 

clear all;
close all;
%% ESPRIT algorithm
% 2019-11-05
% code by wk
%%
dread = pi/180;
N = 10;%阵元个数
M = 3;%信源个数
K = 1024;
snr = 10;
theta =[-30,40,60];%待估角度
lamda = 1;%波长
d = 0:lamda/2:(N-1)*lamda/2;%阵列
A = exp(-1i*2*pi*d.'*sin(dread *theta));
S = randn(M,K);
X = A*S;
X1 = awgn(X,snr,'measured');
%%
Rxx = X1 * X1'/K;
[EV,D] = eig(Rxx);
EVA = diag(D)';
[EVA,I] = sort(EVA);
EV = fliplr(EV(:,I));%按特征值从大到小排序
s = EV(:,1:M);%信号子空间
s1 = s(1:N-1,:);
s2 = s(2:N,:);
wa = s2/s1;%s2/s1 =s2*s1^-1;s2\s1 =s2^-1*s1
[Ew,Dw] = eig(wa);
f= diag(Dw)'; 
source_doa=zeros(1,M);
for j = N-M:N-1
  source_doa(j-(N-M)+1)=asin(angle(f(j))/pi)*180/pi;  
end
disp('source_doa');
disp(source_doa);

%这个算法在看完书之后，理解到就能很快写出代码。

clear all;
close all;
%% TLS-ESPRIT algorithm
% 2019-11-05
% code by wk
%%
dread = pi/180;
N = 10;%阵元个数
M = 3;%信源个数
K = 1024;
snr = 10;
theta =[-50,30,60];%待估角度
lamda = 1;%波长
d = 0:lamda/2:(N-1)*lamda/2;%阵列
A = exp(-1i*2*pi*d.'*sin(dread *theta));
S = randn(M,K);
X = A*S;
X1 = awgn(X,snr,'measured');
Rxx = X1 * X1'/K;
[EV,D] = eig(Rxx);
EVA = diag(D)';
[EVA,I] = sort(EVA);
EV = fliplr(EV(:,I));%按特征值从大到小排序
s = EV(:,1:M);%信号子空间
s1 = s(1:N-1,:);
s2 = s(2:N,:);
s12 = [s1,s2];
[Es,Ds] = eig(s12'* s12);
Ess = diag(Ds)';
[Ess,I] = sort(Ess);
Es = fliplr(Es(:,I));%按特征值从大到小排序,找出最小的特征值约等于0.
C = mat2cell(Es,[M,M],[M M]);
U12 = cell2mat(C(1,2));
U22 = cell2mat(C(2,2));
fi = -1*(U12/U22);
[Ef,Df] = eig(fi);
f = diag(Df)'; 
source_doa=zeros(1,M);
for j = 1:M
  source_doa(j)=asin(angle(f(j))/pi)*180/pi;  
end
disp('doa angle:');
disp(sort(source_doa));
%%这个出现的问题在s12'* s12的特征值分解，只是单纯的分解，但是特征值对应的特征向量有关系，需考虑有M个入射角度，

%但总共2M*2M的矩阵，意味着有M列的噪声，所以按照博客的算法分析，需要将信号矩阵和噪声矩阵排序，找到噪声矩阵。