关于ESPRIT的Matlab例程详解1

学谱估计做ISAR成像，由于是个新坑，一切还得自学，特地把潜水的博客用起来，转一下别人的例子，做下笔记。而且线代矩阵的知识早就还给老师了，一切从头补吧（泪）。

原文链接

ESPRIT算法

ESPRIT算法，estimating signal parameter via rotational invariance techniques，中文名“基于旋转不变技术的信号参数估计”。

同样的，和Musci,Root-Music算法一样，这类基于相关矩阵特征分解的信号频率估计方法，书上都是基于复正弦加白噪声信号模型： 

 

 

 

根据M个时刻的观测值可以表示出向量形式，x(n)=As(n)+v(n)；然后定义y(n)=x(n+1)，同理写出矩阵形式，并找出y(n)和s(n)的关系：

 

 这里的φ是e^(-jwK)组成的对角矩阵。

求x(n)的自相关(与y(n)的自相关是相等的)和与y(n)的互相关；然后对Rxx进行特征值分解，然后定义两个矩阵：

 

 

 对上面那个矩阵对放在一起并求广义特征值，然后由非零解和矩阵是否奇异得出秩的关系:

 

 

所以广义特征值恰为φ^H的对角元，其与信号频率有关。

 

 

                                                            

ESPRIT算法代码：在代码中由于相关矩阵都是估计的，对矩阵对进行广义特征值分解时，取的是最接近圆的K个广义特征值。K是由输入向量里有多少个频率得出的。

clear all,close all
N=1000;%样本数
noise=(rand(1,N)+1i*rand(1,N))/sqrt(2);
n=1:1:N;
un(n)=exp(1i*0.5*pi*n+1i*2*pi*rand(1))+exp(-1i*0.3*pi*n+1i*2*pi*rand(1))+noise;%实际成像是要把接收到的信号化为全极点模型，这个全极点模型所包含的正弦信号的频率就是我们要估计的值，这些频率值将通过模型间的关系与散射点坐标联系起来。
%% 构造自相关矩阵
M=8;
for k=1:N-M
xs(:,k)=un(k+M-1:-1:k).'; %从 k+m-1 一直减到 k 并变为列向量  %另外说一句，这里原信号长度是1000，按代码的意思是仿照书上构造子阵，比如说两个子阵，每个子阵1000-8个阵元（完全是为了仿真数据构造方便），行数按我理解是不同的观测时刻，这个数目不宜过少，
　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　%单次的回波信息理论上只能估计出一个散射点的位置信息，且不一定准确。
end
Rxx=xs(:,1:end-1)*xs(:,1:end-1)'/(N-M-1);
Rxy=xs(:,1:end-1)*xs(:,2:end)'/(N-M-1);
%% 利用奇异分解找到最小特征值
[V,D]=svd(Rxx);              %这里svd只有两个返回值的原因是Rxx是中心对称矩阵，奇异值分解得到的矩阵是V*D*V',即V和U是相同的。
ev=diag(D);emin=ev(end);     %diag是求求矩阵的对角元素，并且从大到小排起来，这里需要找到最小值，即最后一个。
%% 构造矩阵对,并对其进行特征分解
z=[zeros(M-1,1),eye(M-1);0,zeros(1,M-1)];
Cxx=Rxx-emin*eye(M);
Cxy=Rxy-emin*z;                   %对应的数学原理看书
[V,D]=eig(Cxx,Cxy);          %求Cxx Cxy组合起来的矩阵的全部特征值，构成对角阵D，并求A的特征向量构成V的列向量。
z=diag(D);                   %广义特征值恰为φ^H的对角元，其与信号频率有关。
w=angle(z)/(2*pi);           %其实是求的f，2pif=w。
err=abs(abs(z)-1);           %这里是个关键点，前提是假设已知散射点数目为2.若无噪声，极点都分布在单位圆上，但这里加了噪声，就去寻找离单位圆最近的极点，方法就是对模值减1再取绝对值。
[Enew,ad] = sort(err); %对其排序，小到大排列
w(ad([1,3]))                 %12和13的效果差不多。