esprit.txt

该程序利用matlab软件仿真了阵列信号处理中的ESPRIT算法。

%--------------------- 参数定义---------------------
N= 128; % 采样点数
M = 5; % 待估频率分量个数
SNR = 30; % 信噪比
omega = [-pi/3 -pi/6 0 pi/4 pi/3]; % 待估的频率分量
s = [1 10 20 10 1]; % 待估频率分量的振幅

%------------------ 采样数据生成(加高斯白噪声) ---------------------
x = zeros(1,N);
for k = 0:(N-1)
x(k+1) = s*exp(j*omega*k).';
end
y = awgn(x,SNR); % 加性白噪声

%------------------ 计算自相关函数r(0),r(1),...,r(m) -----------------------
r = xcorr(y,'biased');
r = conj(r); % xcorr计算的是E[x(n)*y(n-m)]
r(2*N) = 0; % r(m)=0

%------------------ 构造自相关矩阵Rxx和互相关矩阵Rxy ----------------------
Rxx = zeros(N,N);
Rxy = zeros(N,N);
for k=1:N
Rxx(k,:) = r(N-k+1 : 2*N-k);
Rxy(k,:) = r(N-k+2 : 2*N-k+1);
end

%------------------- 计算Rxx的特征值,提取噪声功率σ^2 ------------------------
d = eig(Rxx);
sigma = min(real(d));

%------------------- 计算Cxx和Cxy ---------------------------
Z = circshift(eye(N),[0,-1]);
Z(1,N) = 0; % first subdiagonal matrix
Cxx = Rxx - sigma*eye(N);
Cxy = Rxy - sigma*Z;

%------------------- 计算Cxx和Cxy的广义特征值------------------
gamma = eig(Cxx,Cxy);

%------------------- 在复平面上比较给定频率和估计频率 ---------------------
plot(gamma,'x') % 估计频率
hold on
plot(exp(j*omega),'ro') % 给定频率
hold off