xinyuanshumugde_page46.m

阵列信号处理的几个程序

%%%%%%信源数目估计----盖而圆半径法page46－48
clc;
clear ;
close all;
tic
M=16;%阵列的天线数
N=4;%信源数
snap=1000;%快拍数目
L=snap;
C=3e8;
lamda=0.2;
f0=C/lamda;
d=0.5*lamda;
% k=d/lamda;
theta0=5;
theta1=20;
theta2=40;
theta3=15;
fs=1000;
ts=1/fs;
t=(0:snap-1)*ts;
a=[0:M-1]';%阵列矢量
u0=5;
u1=10;
u2=20;

%相干信号源 
s0=exp(j*2*pi*(f0*t+0.5*u0*t.^2));
s1=exp(j*2*pi*(f0*t+0.5*u0*t.^2));
s2=exp(j*2*pi*(f0*t+0.5*u0*t.^2));
s3=exp(j*2*pi*f0*t);


%阵列流行矢量

a_theta0=exp(j*2*pi*d/lamda*a*sin(theta0/180*pi));
a_theta1=exp(j*2*pi*d/lamda*a*sin(theta1/180*pi));
a_theta2=exp(j*2*pi*d/lamda*a*sin(theta2/180*pi));
a_theta3=exp(j*2*pi*d/lamda*a*sin(theta3/180*pi));

A=[a_theta0 a_theta1 a_theta2 a_theta3];%子阵的导向矢量，每个子阵完全相同
S=[s0;s1;s2;s3];
X0=A*S;
R_s=S*S';

snr=10;%信噪比

%产生满足一定信躁比的噪声
randn('state',0);
real_noise=randn(size(X0));
randn('state',3);
imag_noise=randn(size(X0));
noise=(real_noise+j*imag_noise)/2^0.5;
noise=10^(-snr/20)*noise;

X=X0+noise;
R_x=X*X'/length(t);
R_pie=R_x(1:M-1,1:M-1);
[V,D]=eig(R_pie);
T=[V zeros(M-1,1);zeros(1,M-1) 1];%%%%%%%%构造矩阵T
Rt=T'*R_x*T;
Rou=abs(Rt(1:M-1,M));
temp0=sum(Rou);

for n=0:M-1    %%%%%%%%%M次迭代再求平均
    for k=1:M-1
        dl=0.5;
        GDE(n+1,k)=Rou(k)-dl/(M-1)*temp0;%%%%%%%%%式（2.5.25）
    end
    K0=find(GDE(n+1,:)>=0);%%%%%%%%%找出第一个大于零的数所对应的位置
    K(1,n+1)=M-K0(1,1);       %%%%%%%%%信源数目的估计
end
figure(1);plot(1:M-1,GDE(1:M,:));grid on;
figure(2);plot(0:M-1,K,'.');grid on;
disp('M次迭代后信源数目估计的平均值');
mean(K) %%%%%%%%%%%%平均的结果显示出来