mvdrzhen.m

波束形成中MVDR波束形成的数据仿真及其结果

clear all;
tic
%=======================================================
%==================阵列信号的MVDR波速形成================
%=======================================================
%根据MVDR波束形成理论，对阵列信号波束形成
%涉及参数有，阵列的阵元个数N，阵元间距d,信噪比SNR，积分时间T，带宽B，
%以及采样频率rate，数据分段数M
%
%2007.12  xiaoleng

M=40;%mvdr 分段数

%设定阵元个数，阵元间距及阵元位置
N=5;
d_lambda=1/2;
R_lambda=[zeros(1,N);d_lambda*(0:N-1);zeros(1,N)];
%设定采样频率及数据长度
rate=3200;
length_data1=1024;
length_data0=512;
length_data=(M-1)*length_data0+length_data1;
t=(0:(length_data-1))/rate;


%设定信号比
SNR=5;
SNR_nb=SNR+10*log10(length_data/2);
An=0.707/sqrt(10^(SNR/10));

%入射波方位
Theta=135*pi/180;
Phe=45*pi/180;

%单矢量水听器阵列流型矢量
u=[cos(Theta)*cos(Phe) sin(Theta)*cos(Phe) sin(Phe)].';
h=[1 u.'].';
%矢量水听器阵列的阵列流型矢量
ap_e=exp(j*2*pi*R_lambda.'*u);
a_e=kron(ap_e,h);

%设定信号
%CW信号
f=400;
ps=cos(i*2*pi*f*t);

% %线性调频信号
% f1=400
% f2=800;
% ps=chirp(t,f1,1,f2);

pn=An*randn(1,length_data);
p=ps+pn;
vs=[ps*u(1,1);ps*u(2,1);ps*u(3,1)];
vn=An/sqrt(3)*randn(3,length_data);

x=[p;vs+vn];                        %单矢量水听器的输出数据
X=kron(ap_e,(fft(x.')).');                     %矢量水听器阵列的输出数据



%分段进行MVDR处理
%窄带分割成M段
for m=1:M
XX(:,:,m)=fft(X(:,(m-1)*length_data0+1:(m-1)*length_data0+length_data1).');
end

%进行MVDR处理

%对于窄带仿真信号
fs1=398;
fs2=402;
% 
% %对于宽带信号
% fs1=398;
% fs2=802;

%取出信号频带内的数据
for i=1:M
    for fs=fs1:fs2;
    
    N_fs=fix(fs*length_data/M/rate)+1;
    
    XXX(:,i,fs-fs1+1)=XX(N_fs,:,i);    %XXX(:,:,fs)为各个频点的数据
    end
end
   

for fs=fs1:fs2
    fs
    R=XXX(:,:,fs-fs1+1)*XXX(:,:,fs-fs+1)';  %对于单频点，生成数据协方差矩阵
    c=inv(R);                               %求逆矩阵
    for i=1:180
        phe_x(i)=(i-90)/180*pi;
        for k=1:360
            
            theta_x(k)=k/180*pi;
            uu=[cos(theta_x(k))*cos(phe_x(i)) sin(theta_x(k))*cos(phe_x(i)) sin(phe_x(i))].';
            hh=[1 uu.'].';
            ap_ee=exp(j*2*pi*R_lambda.'*uu);
            a_ee=kron(ap_ee,h);
            
            P_mvdr_t(fs-fs1+1,i,k)=1/(a_ee'*c*a_ee);     %求出对应频点对应方位的强度值
        end
    end
end

P_mvdr(:,:)=mean(P_mvdr_t(1:end,:,:));                        %对各个频点做平均
[theta,phe]=meshgrid(theta_x,phe_x);
mesh(theta,phe,real(P_mvdr)./max(max(real(P_mvdr))));
xlabel('水平方位角');
ylabel('俯仰角');
zlabel('相对强度');

    
toc