%基本 ESPRIT 算法,第二种方法 最小二乘法 clear all;close all;clc; 即 阵元数 c=3*10^8; f=3*10^9; %% 求得信号的波长 lamda=c/f; %% 阵元的间距 d=lamda/2; %% （n-1）为子阵列的个数 n=10; %% 信号的数目 signal_number=3; %% 三个信号的角度值 thita1=-25; thita2=30; thita3=65; %% 三个信号的中心频率 f1=40; f2=20; f3=70; %% 在时域来说，是快拍数（一段时间内对阵列数据采样的个数）；在频域来说，是 DFT 的时间 子段的个数。 snapshot=1:2000; %% S 是信号空间，有三个信号组成 S1=2.72*exp(j*2*pi*f1*snapshot/length(snapshot)); S2=4.48*exp(j*2*pi*f2*snapshot/length(snapshot)); S3=7.37*exp(j*2*pi*f3*snapshot/length(snapshot)); S=[S1;S2;S3]; %% 子阵 1 A1=exp(-j*2*pi*d*[0:n-1]*sin(thita1*pi/180)/lamda).'; A2=exp(-j*2*pi*d*[0:n-1]*sin(thita2*pi/180)/lamda).'; A3=exp(-j*2*pi*d*[0:n-1]*sin(thita3*pi/180)/lamda).'; A=[A1,A2,A3]; %% 噪声 假设为高斯白噪声，均值为零的 N= wgn(10,2000,3); %% 求信噪比的 S1,S2,S3 信噪比依次是 10 20 30 s_power1=10*log(2.72^2/2); s_power2=10*log(4.48^2/2); s_power3=10*log(7.37^2/2); snr1=s_power1-3; snr2=s_power2-3; snr3=s_power3-3; %% 整个阵列接收到的数据 0-n-1 为阵列 1； 1-n 为阵列 2 的 X=A*S+N; %% 协方差矩阵 Rxx=X*X'/length(snapshot); %% 对整个数据的协方差矩阵进行特征分解，从而得到特征值向量 D 和特征向量 V 

[V,D]=eig(Rxx); %[Y,I]=sort(diag(D)); Us=V(:,n-signal_number+1:n); %% 两个方阵张成的两个子空间 U1=Us(1:n-1,:); U2=Us(2:n,:); %% 利用最小二乘法 求得旋转不变关系矩阵，然后进行特征分解 [p,q]=eig(inv(U1'*U1)*U1'*U2); %% 利用上面求得的矩阵来获得角度 for i=1:signal_number; alpha(i)=real(asin(-j*(log(q(i,i)))*lamda/(-2*pi*d))*180/pi); %张贤达《矩阵分析与应用》 第 528 页 end; %% 作图 stem(alpha,ones(1,signal_number),'r--');grid; axis([-90 90 0 2]); text(alpha(1)-4,1.1,num2str(alpha(1)));text(alpha(1)-15,1.4,'信号 1，信噪比为 10'); text(alpha(2)-4,1.1,num2str(alpha(2)));text(alpha(2)-15,1.4,'信号 2，信噪比为 20'); text(alpha(3)-4,1.1,num2str(alpha(3)));text(alpha(3)-15,1.4,'信号 3，信噪比为 30'); ylabel('DOA 估计的角度值 '); xlabel('角度'); title('ESPRIT 算法 DOA 估计');