ESPRIT算法(最小二乘法)matlab程序

%基本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); %张贤达《矩阵分析与应用》第528页%% 利用上面求得的矩阵来获得角度

for i=1:signal_number;

alpha(i)=real(asin(-j*(log(q(i,i)))*lamda/(-2*pi*d))*180/pi);

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估计');