music 方位估计 实验报告三 MATLAB 代码

实验报告三

实验目的：

实现常规波束形成及基于MUSIC 方法的方位估计。 实验内容：

1）若干阵元的接收阵，信号频率为10KHz ，波束主轴12度，仿真给出常规波束形成的波束图。 2）16个阵元的均匀线列阵，信号频率为10KHz ，信号方位为12度，用MUSIC 方法完成目标定向，信噪比-5dB ，0dB ，5dB 。

i) 波束形成时的阵型设计为两种，一种是均匀线列阵，阵元16个；一种是均匀圆阵，阵元数为16个，比较这两种阵型的波束图。

ii ）比较不同信噪比下MUSIC 方法估计的性能（统计100次）。 实验原理：

i)常规波束形成：

如图所示，基阵的输出),(θt v 。

∑∑=*=*

==M

m i i M m i i

w t x t x w t v 1

1

)

()()()(),(θθθ

采用向量符号则有，

)()()()(),(H H θθθw x x w t t t v ==

式中，x(t)和w(q )分别为观测数据向量和加权系数向量，

)

,(θt v 图 1 波束形成器基本原理图

T M 21])()()([)(t x t x t x t =x T M 21])()()([)(θθθθw w w =w

基阵输出端的空间功率谱表示为：

)

()( )()]()([)( )]()()()([ )],(),([ ]),([)(H H H H H *2

θθθθθθθθθθRw w w x x w w x x w =====t t E t t E t v t v E t v E P

式中，R 为观测数据的协方差矩阵。

ii ）基于MUSIC 方法的方位估计：

)()()()(1

t n t s a t x i d

i +=∑=θ

T M 21])()()([)(t x t x t x t =x

)()()()(t n t s A t x +=θ

假设:

(1 ) 信号源的数目d 是已知的, 且d < M ;

(2 ) 各信号的方向矢量是相互独立的, 即)(θA 是一个列满秩矩阵;

(3 ) 噪声)(t n 是空间平稳随机过程, 为具有各态历经性的均值为零、方差为σ2n 的高斯过程;

(4 ) 噪声各取样间是统计独立的。

在上述假设条件下, 基阵输出的协方差矩阵可表示为：

I A AR t x t x E R H s H 2])()([α+==

其中, R s 为信号的协方差矩阵；I 为单位矩阵。对R 进行特征分解, 并以特

征值降值排列可得

H

m

m M

d m m H

m m d

m m e

e e e R ∑∑+==+

=1

1λ

λ

信号子空间与噪声子空间正交。 若噪声子空间记为E N , 即

∑+==

M

d m H

m

m N e

e

E 1

则有：

0)(=i N a E θ

因此，对矩阵R 进行特征分解后, 取噪声特征矢量, 就可得到信号源的到达角。

定义多重信号分类法的空间谱函数为：

)

()(1

θθa E E a P H N

N H

MUSIC =

实验过程及结论：

1.常规波束形成，阵列为均匀线列阵：

-100

-80-60-40-20

020*********

-70-60

-50

-40

-30

-20

-10

均匀线列阵 12度主轴

方位角/度

归一化功率/d B

2.常规波束形成，阵列为均匀圆阵：

-200

-150-100-50

050100150200

-80-70-60-50-40-30-20-10

0均匀圆阵 12度主轴

方位角/度

归一化功率/d B

结论：由上面两个仿真结果可以看出，采用均匀线列阵和均匀圆阵都可以完成波束形成，都可以获得阵处理增益。

3.MUSIC

算法方位估计

：

-100

-80-60-40

-20020406080100

-60-50

-40

-30

-20

-10

方位角 / 度

归一化功率 / d B

结论：从上图得出，信噪比越高，MUSIC方位估计越准确。

附录：MATLAB仿真程序代码1.常规波束形成

%%

clear all;close all;clc

%%%%%%常规波束形成

%信号参数

f0=10e3;

fs=10*f0;

T=2e-3;

t=1:1/fs:T-1/fs;

%均匀线列阵阵列参数

c=1500;

lamda=c/f0;

d=lamda/2;

M=16;

direct=12;%波束指向

%阵列流行

theta=-90:90;

for i_theta=1:length(theta)

fai=2*pi*d*sin(theta(i_theta)*pi/180)/lamda;

for i_M=1:M

A(i_M,i_theta)=exp(-sqrt(-1)*(i_M-1)*fai);

end

end

%加权向量w

fai_direct=2*pi*f0*d*sin(direct*pi/180)/c;

for i_M=1:M

w(i_M)=exp(sqrt(-1)*(i_M-1)*fai_direct);

end

w=w';

%无噪声输出

P=w'*A;

P=20*log10(abs(P)/max(abs(P)));

figure(1)

plot(theta,P);