ESPRIT算法记录

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

前言

提示：这里可以添加本文要记录的大概内容：

ESPRIT算法记录

---

提示：以下是本篇文章正文内容，下面案例可供参考

一、ESPRIT算法什么？

ESPRIT算法是一种基于子空间的角度估计算法
对于均匀线阵模型，显然有
$$Y=Ax+n$$
其中A是导向矢量构成的矩阵，$A=[a({\theta }_1),a({\theta }_2)...,a({\theta }_k)]$,$a(\theta )=[e^{i\pi \ast 0\ast sin(\theta )},..,e^{i\pi \ast (N-1)sin(\theta )}{]}^T$,其中N为阵元数，k为用户数。

二、ESPRIT算法

1.计算理论协方差矩阵

$C=E[Y{Y}^H]=AP{A}^H+{\sigma }^2$，P为一个对角阵，所以可以得到A构成了信号子空间。

2.旋转不变性

可以看出A具有一些特性，即A的每一行和上面都是成比例的，所以我们有
$$J_{1}Adiag(e^{i\pi \ast \Delta \ast sin({\theta }_l)}{)}_{l-1}^k=J_2\ast A$$,其中$J_1,J_2$定义为从N行中取n行，$J2$与$J1$相差行数为$\Delta$，实际上主要是A导向矢量的特性。
实际上如果我们能够估计的很好那么显然$Us=A\ast Q$，因为如果估计的很好$Us$和$A$就是一个空间，所以存在一个变换，将两者联立起来，Q是一个满秩矩阵。

所以我们有下面式子
$$Qdiag(e^{i\pi \ast \Delta \ast sin({\theta }_l)}{)}_{l-1}^k{Q}^{-1}=(J_1{U}_s{)}^{+}{J}_2{U}_s$$

其实也就是$diag(e^{i\pi \ast \Delta \ast sin({\theta }_l)}{)}_{l-1}^k$与$(J_1{U}_s{)}^{+}{J}_2{U}_s$，相似，而我们要求的角度就在左边这个矩阵的对角线元素上（其实也就是特征值），而由于相似，我们可以直接去求右边矩阵的特征值，然后通过反变换
$$\theta =arcsin(arg(\lambda ((J_1{U}_s{)}^{+}{J}_2{U}_s))/(pi\ast \Delta )),$$

---

3.Code

%% ESPRIT theory, k =1 error

close all; clear; clc

coeff = 1;
N = 100*coeff;
T = 2*N*coeff;
c= N/T;
p = N;
e = @(index) [zeros(index-1,1);1;zeros(p-index,1)];
% theta_true = -10/180*pi; % should belong to (-pi, pi)

clear i
theta_true = [10];
k = length(theta_true);
A = exp(pi*1i*(0:N-1)'*sind(theta_true))/sqrt(N);
P = diag([1]);
S = 1*randn(k,T);

Z = complex(randn(N,T), randn(N,T))/sqrt(2); %% 噪声功率为1
X = A*S + Z;

SCM = X*(X')/T;

[U,eigs_SCM] = eig(SCM,'vector');
[eigs_SCM, index] = sort(eigs_SCM,'descend');
U = U(:, index);
u_S = U(:,1:k);

J_tmp = eye(N);
n = round(N/2);
Delta = 5;
index = 30;

J1 = J_tmp(index:n+index,:);
J2 = J_tmp(index+Delta:n+index+Delta,:);

test1 = u_S'*(J1'*J2)*u_S;
test2 = u_S'*(J1'*J1)*u_S;


test1_estim = A'*(J1'*J2)*A;
test2_estim = A'*(J1'*J1)*A;

Res = inv(test2_estim) * test1_estim;
[EigenVector,EigenValue] = eig(Res,'vector');
res1 = asind(angle(EigenValue)/pi/Delta)
figure;
plot(res1,0,'.','Color','g','MarkerSize',30)

hold on;
Res = inv(test2) * test1;
[EigenVector,EigenValue] = eig(Res,'vector');
res2 = asind(angle(EigenValue)/pi/Delta)
plot(res2,0,'.','Color','r','MarkerSize',30)
% xline(theta_true,0,'-')
% plot(theta_true,0,'.','Color','b','MarkerSize',30)

给出了协方差矩阵估计正确和误差近似的结果，因为实际上我们只能得到取样协方差矩阵，而不是真正的协方差矩阵，所以我们用${\tilde{U}}_s$去替代真实的$U_s$，是有误差的，图上给出了误差实例

多角度

%% Estimation of the "direction of arrival" of signal from noisy observations 
% MUSIC versus ESPRIT
close all; clear; clc

coeff = 2;
p = 256*coeff;
n  = 3*p*coeff;

theta_true = [-10, 35, 37]./180*pi; % should belong to (-pi, pi)
%theta_true = [35, 37]./180*pi; % should belong to (-pi, pi)
%theta_true = [-10, 30]./180*pi; % should belong to (-pi, pi)
k = length(theta_true);
sigma2 = .1;
%P = eye(k);
P = diag([1 3 7]);

clear i

a = @(theta) exp(pi*1i*sin(theta)*(0:p-1)')/sqrt(p);
%a = @(theta) exp(pi*1i*sin(theta)*(0:p-1)');

A = [];
for tmp_index=1:length(theta_true)
    A = [A a(theta_true(tmp_index))];
end

%theta_range = linspace(-180,180,500)./180*pi;
theta_range = linspace(-90,90,1000)./180*pi;


S = sqrtm(P)*randn(k,n);
Z = complex(randn(p,n), randn(p,n));
X = A*S + sqrt(sigma2/2)*Z;

SCM = X*(X')/n;

[U,eigs_SCM] = eig(SCM,'vector');
[eigs_SCM, index] = sort(eigs_SCM,'descend');
U = U(:, index);
U_S = U(:,1:k);

J_tmp = eye(p);
n = round(p/8);
J1 = J_tmp(1:n,:);
J2 = J_tmp(2:n+1,:);

test1 = U_S'*(J1'*J2)*U_S;
test2 = U_S'*(J1'*J1)*U_S;
Psi_test = inv(test2) * test1;
[EigenVector,EigenValue] = eig(Psi_test,'vector');
resU = asind(angle(EigenValue)/pi/1)

test1 = A'*(J1'*J2)*A;
test2 = A'*(J1'*J1)*A;
Psi_test = inv(test2) * test1;
[EigenVector,EigenValue] = eig(Psi_test,'vector');
resA = asind(angle(EigenValue)/pi/1)

总结

记录学习过程