chapter5_10.m

阵列信号处理中各中方位来波方向估计

%%%%chapter5_10%%%%%
%%%利用投影矩阵法进行空间谱估计
% %信噪比5db感觉太小了,得出的结果不够理想
clc;clear;
P=2;%信号源个数
fai=[30,40]*pi/180;%信号的方向转为弧度
N=100;%N为快拍得数目
fai=pi*sin(fai);
theta=linspace(-100,100,1000);%画图时的角度范围,并有1000个点
L=8;%为阵元数目

%%%产生A_theta%%%%
A_theta=zeros(L,2);
for k=1:L
    A_theta(k,1)=exp(j*fai(1)*(k-1));
    A_theta(k,2)=exp(j*fai(2)*(k-1));
end

%%%% 产生信号%%%%
Sig=zeros(2,N);
for k=1:N
    Sig(1,k)=exp(j*pi/4*(k-1));
    Sig(2,k)=exp(j*pi/6*(k-1));
end
SNR=[-5,0,5];%分别是三个信噪比
R_guji=zeros(L,L);
y_gce=zeros(L,1);%实际的观测值
sums=zeros(L,L);
for k=1:N
    y_gce=1/2^0.5*sqrt(10.^(SNR(3)/10))*A_theta*Sig(:,k)+randn(L,1);%SNR(i)分别代表三个信噪比的取值
    sums=sums+y_gce*y_gce';
end
R_guji=sums/N;
%%%%利用投影矩阵法进行估计%%%%%

%%%(1)取得是自相关矩阵的第一,二列作为投影矩阵%%%%
I=eye(8,8);
Q1=zeros(8,2);
P_pr1=zeros(1,1000);
P_pr2=zeros(1,1000);
Q1=R_guji(:,1:2);%Q1为构造的投影矩阵
P_1=I-Q1*inv(Q1'*Q1)*Q1';%求出的是阵列的噪声子空间的投影矩阵;

%%%%(2)取得是自相关矩阵的第一,和最后一列;
Q2=zeros(8,2);
Q2=R_guji(:,[1,end]);
P_2=I-Q2*inv(Q2'*Q2)*Q2';%阵列的噪声子空间的投影矩阵

%%%%谱估计的图形绘制%%%
 theta=linspace(-100,100,1000);  
 theta1=theta*pi/180;
 a_theta=zeros(L,1);
 for k=1:1000
     for n=1:L
         a_theta(n,1)=exp(j*pi*sin(theta1(k))*(n-1));
     end
     P_pr1(k)=1/((a_theta)'*P_1*a_theta);
     P_pr2(k)=1/((a_theta)'*P_2*a_theta);
 end
figure
plot(theta,abs(P_pr1))
xlabel('角度');
ylabel('幅度');
title('投影矩阵法,以第一,二列为噪声子空间')
grid on 
figure
plot(theta,abs(P_pr2))
xlabel('角度');
ylabel('幅度');
title('投影矩阵法,以第一和最后一列为噪声子空间')
grid on 
%%%(3)取第一和第二列,并避开对角线元素%%%

%%%R_guji为8x8的矩阵,所以可取前两列的第三行到第八行%%%%
%%%这时的a_theta为L-2维的列向量%%%
I_2=eye(L-2);
Q3=R_guji(P+1:L,1:2);
P_3=I_2-Q3*inv(Q3'*Q3)*Q3';%阵列的噪声子空间的投影矩阵
%%%%绘制图形%%%%%
a_theta1=zeros(L-2,1);
 for k=1:1000
     for n=1:L-2
         a_theta1(n,1)=exp(j*pi*sin(theta1(k))*(n-1));
     end
     P_pr3(k)=1/((a_theta1)'*P_3*a_theta1);%此为谱估计
 end
figure
plot(theta,abs(P_pr3))
grid on 
xlabel('角度');
ylabel('幅度');
title('投影矩阵法,以第一,二列矢量并避开对角线元素')