最大信噪比准则MATLAB实现

clear all
close all
clc
N=8;%阵元个数
c=3e+8;%c为光速
w=500e6;%雷达工作频率
wavelength=c/w;%波长
%*************以下内容为可修改参数*************%
power0=1;%信号功率
power1=10;%噪声功率
da=wavelength/2;%阵元间距
%da=3;
theta=-(90/180)*pi:0.005:(90/180)*pi;
theta_0=0;%来波方向
theta_1=5/8*pi; %干扰方向
%theta_1=1;
%设置导向矢量,将导向矢量放在一个矩阵中，每一列为一个接受阵元的导向矢量
for i=1:N
a(:,i)=exp(j*2*pi*(i-1)*sin(theta)*da/wavelength);
a_theta0(i)=exp(j*2*pi*(i-1)*sin(theta_0)*da/wavelength);%此处是行向量
a_noise(i)=exp(j*2*pi*(i-1)*sin(theta_1)*da/wavelength);
end
%以下是设置导向矢量的另一种方法
% for i=1:N
% for k=1:length(theta)
% a(k,i)=exp(j*2*pi*(i-1)*sin(theta(k))*da/wavelength);
% end
% end
Rs=power0*a_theta0.'*conj(a_theta0);
Rn=power1*a_noise.'*conj(a_noise)+(5e-2)*eye(N,N);%干扰加噪声，消除干扰于信号之间的相干性
%W=Rs*pinv(Rn)*a_theta0';%pinv(Rn)是矩阵求逆
%以下实现寻找矩阵Rs*pinv(Rn)的最大特征值及其对应的特征向量
%具体方法是：找出最大的特征值，然后用特征值组成的向量减去该值，此时最大特征值处为零，
%让后对得到的向量取非，只有原最大特征值处为一，其余为零，再用find函数得到其位置
% [V,D]=eig(Rn\Rs);
[V,D]=eig(pinv(Rn)*Rs);
[maxr l]=max(diag(D));
% vv=max(D)-lanmda;
% [i,j,z]=find(~abs(vv));%得到的j就是最大特征值及其对应的特征向量的坐标
W=V(:,l);
%W=(Rs./Rn)*a_theta0';
dir=W'*a.';
%theta=-pi/2:0.05:pi/2;
figure(1)
plot(theta*180/pi,abs(dir/max(dir)))
title('SNR方向图')
grid on;
Rx=Rs+Rn;
%实现有E[y(t).^2]=W'*Rx*W,但是此处的a为行向量，W是列向量，所以表达方式不同
ener=conj(a)*Rx*a.';
energy=sum(ener);
% ener=zeros(1,length(theta));
% for i=1:length(theta)
% ener=ener+energy(:,i);
% end
% figure(2)
% plot(theta,20*log(energy/max(energy)))
% title('SNR功率谱')
% grid on