Capon 波束形成matlab仿真(附源代码)

合集下载

CAPON波束形成_Matlab程序

C A P O N波束形成_M a t l a b程序(总3页)--本页仅作为文档封面，使用时请直接删除即可----内页可以根据需求调整合适字体及大小--CAPON波束形成器仿真1．实验原理波束形成就是从传感器阵列重构源信号。

（1）、通过增加期望信源的贡献来实现；（2）、通过抑制掉干扰源来实现。

经典的波束形成需要观测方向(期望信源的方向)的知识。

盲波束形成试图在没有期望信源方向信息的情况下进行信源的恢复。

波束形成技术的基本思想是：通过将各阵元输出进行加权求和，在一时间内将天线阵列波束“导向”到一个方向上，对期望信号得到最大输出功率的导向位置即给出DOA估计。

虽然阵列天线的方向图是全方向的，但阵列的输出经过加权求和后，却可以被调整到阵列接收的方向增益聚集在一个方向上，相当于形成了一个”波束”。

这就是波束形成的物理意义所在。

在智能天线中，波束形成是关键技术之一，是提高信噪比、增加用户容量的保证，能够成倍地提高通信系统的容量，有效地抑制各种干扰，并改善通信质量。

波束形成器的最佳权向量w取决于阵列方向向量)(ak，而在移动通信里用户的方向向量一般未知，需要估计(称之为DOA估计)。

因此，在计算波束形成的最佳权向量之前，必须在已知阵列几何结构的前提下先估计期望信号的波达方向。

Capon 波束形成器求解的优化问题可表述为w arg min P(w)θ=其约束条件为 1)(a w H =θCapon 波束形成器在使噪声和干扰所贡献的功率为最小的同时，保持了期望信号的功率不变。

因此，它可以看作是一个尖锐的空间带通滤波器。

最优加权向量w 可以利用Lagrange 乘子法求解，其结果为 )(a R ˆ)(a )(a R ˆw 1H 1CAP θθθ--= 当μ不取常数，而取作)(a R ˆ)(a 11H θθμ-=时，最佳权向量就转变成Capon 波束形成器的权向量。

空间谱为)(a R ˆ)(a 1)(P 1-H CAP θθθ=2．变量定义 M ：均匀线阵列数目P ：信号源个数nn ：快拍数angle1、angle2、angle3：信号来波角度u ：复高斯噪声Ps ：信号能量refp ：信噪比（实值）X：接收信号Rxx：接收信号的相关矩阵doa：波达方向估计3．仿真结果采用上述算法进行仿真，结果如图所示。

波束形成Matlab程序

1?均匀线阵方向图%8阵元均匀线阵方向图，来波方向为clc; clear all; close all;0度imag=sqrt(_1);element_num=8;% 阵元数为8d_lamda=1/2;%阵元间距d与波长lamda的关系theta=li nspace(-pi/2,pi/2,200);theta0=0;% 来波方向w=exp(imag*2*pi*dl_lamda*sin(theta0)*[0:eleme nt_nu m-1]');for j=1:le ngth(theta)a=exp(imag*2*pi*dd_l amda*si n(theta(j))*[0:eleme nt_nu m-1]'); p(j)=w'*a;end figure;plot(theta,abs(p)),grid on xlabel('theta/radia n')ylabel('amplitude')title('8 阵元均匀线阵方向图')°2 8阵元均匀线阵方向图7654321-15 -1 -0 5 0 06thetaradian1 158当来波方向为45度时，仿真图如下8阵元均匀线阵方向图如下，来波方向为0 度，20log (dB)8阵元均苛銭阵方向图来波方向为0度随着阵元数的增加，波束宽度变窄，分辨力提高：仿真图如下Q d pE =ro 二2. 波束宽度与波达方向及阵元数的关系clcclear allclose allima=sqrt(-1);element_num1=16; %阵元数element_num2=128;element_num3=1024;lamda=0.03; %波长为0.03 米d=1/2*lamda; % 阵元间距与波长的关系theta=0:0.5:90; forj=1:length(theta);fai(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_n um1*d));psi(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_n um2*d));beta(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_ num3*d));endfigure; plot(theta,fai,'r',theta,psi,'b',theta,beta,'g'),grid on xlabel('theta'); ylabel('Width in radians')title(' 波束宽度与波达方向及阵元数的关系') 仿真图如下：。

波束赋形仿真代码

以下是一个简单的波束赋形仿真代码，使用Python编写：python复制代码import numpy as npimport matplotlib.pyplot as plt# 设置参数num_sensors = 10# 传感器数量num_signals = 5# 信号数量num_time_steps = 100# 时间步数# 生成信号signals = np.random.randn(num_signals, num_time_steps)# 生成传感器位置sensor_positions = np.random.rand(num_sensors, 2) * 10# 在10x10的区域内随机生成位置# 生成波束赋形器参数weights = np.random.randn(num_sensors, num_signals)# 进行波束赋形output = np.dot(sensor_positions, weights.T)# 绘制结果plt.figure()plt.subplot(2, 1, 1)plt.imshow(sensor_positions[:, 0], aspect='auto') # 绘制传感器位置图plt.title('Sensor Positions')plt.xlabel('x')plt.ylabel('y')plt.subplot(2, 1, 2)plt.plot(output) # 绘制波束赋形结果图plt.title('Beamforming Output')plt.xlabel('Time Step')plt.ylabel('Amplitude')plt.show()这个代码使用随机生成信号和传感器位置，并使用随机生成的权重进行波束赋形。

输出结果是一个二维图像，其中第一个子图显示了传感器位置，第二个子图显示了波束赋形结果。

自适应波束形成与Matlab程序代码注解

gridon;
xlabel('theta/radian')
ylabel('amplitude/dB')
title([num2str(element_num)'阵元均匀线阵方向图','来波方向为'num2str(theta0*180/pi)'度']);
axis([-1.5 1.5 -50 0]);
(2)仿真结果
A.来波方向为0°
不归一化
归一化
B.来波方向为45°
不归一化
归一化
C.随着阵元数的增加，波束宽度变窄，分辨力提高，仿真图如下：
非归一化
归一化
不归一化
归一化
2.波束宽度与波达方向及阵元数的关系
（1）matlab 程序
clc;
clearall;
closeall;
imag=sqrt(-1);
element_num1=16;
psi(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lambda/(element_num2*d));
beta(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lambda/(element_num3*d));
end
gridon;
xlabel('theta/radian')
ylabel('amplitude/dB')
title([num2str(element_num)'阵元均匀线阵方向图','来波方向为'num2str(theta0*180/pi)'度']);

Matlab波束形成程序

波束形成与智能天线1.均匀线阵方向图%8阵元均匀线阵方向图，来波方向为0度clc;clear all;close all;imag=sqrt(-1);element_num=8;%阵元数为8d_lamda=1/2;%阵元间距d与波长lamda的关系theta=linspace(-pi/2,pi/2,200);theta0=0;%来波方向w=exp(imag*2*pi*d_lamda*sin(theta0)*[0:element_num-1]');for j=1:length(theta)a=exp(imag*2*pi*d_lamda*sin(theta(j))*[0:element_num-1]');p(j)=w'*a;endfigure;plot(theta,abs(p)),grid onxlabel('theta/radian')ylabel('amplitude')title('8阵元均匀线阵方向图')当来波方向为45度时，仿真图如下：8阵元均匀线阵方向图如下，来波方向为0度，20log（dB）随着阵元数的增加，波束宽度变窄，分辨力提高：仿真图如下：2.波束宽度与波达方向及阵元数的关系clcclear allclose allima=sqrt(-1);element_num1=16; %阵元数element_num2=128;element_num3=1024;lamda=0.03; %波长为0.03米d=1/2*lamda; %阵元间距与波长的关系theta=0:0.5:90;for j=1:length(theta);fai(j)=theta(j)*pi/-asin(sin(theta(j)*pi/)-lamda/(element_num1*d));psi(j)=theta(j)*pi/-asin(sin(theta(j)*pi/)-lamda/(element_num2*d));beta(j)=theta(j)*pi/-asin(sin(theta(j)*pi/)-lamda/(element_num3*d)); endfigure;plot(theta,fai,'r',theta,psi,'b',theta,beta,'g'),grid onxlabel('theta');ylabel('Width in radians')title('波束宽度与波达方向及阵元数的关系')仿真图如下：3. 当阵元间距时，会出现栅瓣，导致空间模糊。

波束形成 Matlab程序

1.均匀线阵方向图%8阵元均匀线阵方向图，来波方向为0度clc;clear all;close all;imag=sqrt(-1);element_num=8;%阵元数为8d_lamda=1/2;%阵元间距d与波长lamda的关系theta=linspace(-pi/2,pi/2,200);theta0=0;%来波方向w=exp(imag*2*pi*d_lamda*sin(theta0)*[0:element_num-1]');for j=1:length(theta)a=exp(imag*2*pi*d_lamda*sin(theta(j))*[0:element_num-1]');p(j)=w'*a;endfigure;plot(theta,abs(p)),grid onxlabel('theta/radian')ylabel('amplitude')title('8阵元均匀线阵方向图')当来波方向为45度时，仿真图如下：8阵元均匀线阵方向图如下，来波方向为0度，20log（dB）随着阵元数的增加，波束宽度变窄，分辨力提高：仿真图如下：2.波束宽度与波达方向及阵元数的关系clcclear allclose allima=sqrt(-1);element_num1=16; %阵元数element_num2=128;element_num3=1024;lamda=0.03; %波长为0.03米d=1/2*lamda; %阵元间距与波长的关系theta=0:0.5:90;for j=1:length(theta);fai(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_num1*d));psi(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_num2*d));beta(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_num3*d)); endfigure;plot(theta,fai,'r',theta,psi,'b',theta,beta,'g'),grid onxlabel('theta');ylabel('Width in radians')title('波束宽度与波达方向及阵元数的关系')仿真图如下：3. 当阵元间距时，会出现栅瓣，导致空间模糊。

Capon波束形成器的仿真

Capon 波束形成器的仿真一、原理1、波束形成的定义波束形成就是从传感器阵列重构源信号，然后通过增加期望信源的贡献和抑制掉干扰源来实现的。

实际上波束形成可以看成空域（或时域）滤波器。

2、波束形成器的最佳权向量的推导首先，假定一天线阵列中个阵元的接收信号向量为()n x ，权向量为12[,,,]M w w w =w ，由于阵列的输出是对阵元的接收信号向量在个阵元上个分量的加权和，所以阵列的输出可写作：*1()()()MHmm m n n w x n ===∑y w x再令空间远场期望信号为()d t ，J 个干扰信号为()j i t ，1,2,,j J =，每个阵元上的加性白噪声为()k n t ，他们都具有相同的方差2n σ。

在这些假设条件下，第k 个阵元上的接收信号可以表示为1()()()()()()j Jk k d k i j k j x t a d t a i t n t θθ==++∑111212()()()()()()(),(),()()()()J d i i J M M x t d t n t x t i t n t a a a i t n t x t θθθ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎡⎤=+⎣⎦⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ 矩阵形式如下：1()()()()()()()()j Jd i j j t t t a d t a i t n t θθ==+=++∑x As n由以上可知：N 个快拍的波束形成器的输出为()()(1,2,,)H y t w x t t N ==输出的平均功率为2211222222111111()()()111()()()()()j N N H t t N J NNH Hd j i t j t t P w y t w x t N N w a d t i t w a w n t N N N θθ========⎡⎤=++⎢⎥⎣⎦∑∑∑∑∑∑当N →∞时，上式变为{}{}{}{}22222221()()()()()()()()jHHHJHHdj i nj P w E y t w E x t x t w wRwE d t w a E i t w a wθθσ=====++∑为了保证来自方向d θ的期望信号的正确接收，并完全抑制其它J 个干扰，关于权向量的约束条件应为：()1H d w a θ=和()0j H i w a θ=加上约束条件后{}222()()n P w E d t w σ=+波束形成器最佳权向量的确定以叙述为：在置零约束条件的约束下，求满足{}{}2min ()min HwwE y t w Rw =的权向量w 。

波束形成matlab程序

psi(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_num2*d));
beta(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_num3*d));
end
a=exp(imag*2*pi*d_lamda*sin(theta(j))*[0:element_num-1]'); p(j)=w'*a; end figure; plot(theta,abs(p)),grid on xlabel('theta/radian') ylabel('amplitude') title('8 阵元均匀线阵方向图 ')
当来波方向为 45 度时，仿真图如下：
8 阵元均匀线阵方向图如下，来波方向为 0 度， 20log（ dB ）
随着阵元数的增加，波束宽度变窄，分辨力提高：仿真图如下：
2. 波束宽度与波达方向及阵元数的关系
clc
clear all
close all
ima=sqrt(-1);
element_num1=16;
1.均匀线阵方向图
%8 阵元均匀线阵方向图，来波方向为 0 度 clc; clear all; close all; imag=sqrt(-1); element_num=8;% 阵元数为 8 d_lamda=1/2;% 阵元间距 d 与波长 lamda 的关系 theta=linspace(-pi/2,pi/2,200); theta0=0;% 来波方向 w=exp(imag*2*pi*d_lamda*sin(theta0)*[0:element_num-1]'); for j=1:length(theta)

（完整版）波束形成-Matlab程序

（完整版）波束形成-Matlab程序1.均匀线阵方向图%8阵元均匀线阵方向图，来波方向为0度clc;clear all;close all;imag=sqrt(-1);element_num=8;%阵元数为8d_lamda=1/2;%阵元间距d与波长lamda的关系theta=linspace(-pi/2,pi/2,200);theta0=45/180*pi;%来波方向 (我觉得应该是天线阵的指向)%theta0=0;%来波方向w=exp(imag*2*pi*d_lamda*sin(theta0)*[0:element_num-1]');for j=1:length(theta) %(我认为是入射角度，即来波方向,计算阵列流形矩阵A) a=exp(imag*2*pi*d_lamda*sin(theta(j))*[0:element_num-1]');p(j)=w'*a; %(matlab中的'默认为共轭转置，如果要计算转置为w.'*a)endfigure;plot(theta,abs(p)),grid onxlabel('theta/radian')ylabel('amplitude')title('8阵元均匀线阵方向图')见张小飞的书《阵列信号处理的理论和应用2.3.4节阵列的方向图》当来波方向为45度时，仿真图如下：8阵元均匀线阵方向图如下，来波方向为0度，20log（dB）随着阵元数的增加，波束宽度变窄，分辨力提高：仿真图如下：2.波束宽度与波达方向及阵元数的关系clcclear allclose allima=sqrt(-1);element_num1=16; %阵元数element_num2=128;element_num3=1024;lamda=; %波长为0.03米d=1/2*lamda; %阵元间距与波长的关系theta=0::90;for j=1:length(theta);fai(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_num1*d)); psi(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_num2*d));beta(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_num3*d)); endfigure;plot(theta,fai,'r',theta,psi,'b',theta,beta,'g'),grid onxlabel('theta');ylabel('Width in radians')title('波束宽度与波达方向及阵元数的关系')仿真图如下：3.当阵元间距/2dλ>时，会出现栅瓣，导致空间模糊。

Matlab仿真相控阵波束扫描的“动图”（附代码）

Matlab仿真相控阵波束扫描的“动图”（附代码）在《无源/有源相控阵雷达、电扫阵列及典型系统》中有下面这样一幅动图，有人问是怎么画出来的。

今天我们就来分析一下，并且手把手教大家动手画出来。

动图分析：我们观察到这是一个极坐标下线阵的波束扫描图，其中外圈标注的是度数0~360°，0°表示垂直于线阵的方向。

其中，一圈一圈的环代表的是幅度，标出的100,80,60...是dB值。

波束增益最大的是指向0°的方向，最大值100dB，旁瓣比主瓣低接近40dB，说明进行了降低旁瓣的加权处理，并且波束在扫描到30多度后出现了删瓣。

动图仿真：这里我们做了一定的简化，首先，没有进行降低旁瓣的加权。

幅度进行了归一化，也就是指向0°时是0dB。

仿真时设置：N=32，d=λ/2，波束指向θB从0~60°，再从60~0°，得到下面的动图：从该动图可以看出：1. 随着波束指向的增大，波束宽度变胖，当达到60°时，变胖了一倍。

2. 随着波束指向的增大，增益也有降低。

3. 没有进行降低旁瓣的加权，旁瓣较高。

几个知识点科普：波束宽度与扫描角θB的关系：当扫描的最大角度为θmax时，为了不出现删瓣，阵元间距d和波长λ需要满足关系：也就是说当阵元间距小于半波长时，即使扫描到90°都不会出现删瓣。

如果你想看看出现删瓣的情况，这里也有(设置d=0.7λ)。

有一点需要注意，这里仿真图中90~270°内出现波束是因为这里仿真的每个阵元都是全向天线。

仿真程序如下：Matlab中有函数polarplot可以直接使用，但不是dB表示的。

需要polarplotdb函数以及上面仿真程序的可以发送“180724”直接下载。

常规波束形成matlab

常规的波束形成方法通常使用MATLAB实现。

以下是一个简单的MATLAB代码示例，演示了如何使用常规的波束形成方法：```matlab设置参数c = 343; 光速，单位为m/sf = 1000; 信号频率，单位为HzN = 128; 阵元数目d = c/f; 阵元间距，单位为mtheta = 1:1:180; 方向角，单位为度t = 0:1/f:1-1/f; 时间，单位为sA = zeros(1,length(theta)); 信号幅度A(5) = 5; 在第5个方向上有目标W = zeros(N,length(theta)); 阵列流形计算阵列流形for i = 1:1:Nfor j = 1:length(theta)W(i,j) = exp(-2*pi*f*(i-1)*d*sind(theta(j))/c);endend计算波束形成输出S = zeros(length(theta),length(t));for j = 1:length(theta)S(j,:) = sin(2*pi*f.*t).*A(j);endS = W*S;画图figure;plot(theta,abs(S));xlabel('方向角（度）');ylabel('幅度');title('波束形成输出');```这段代码使用了常规的波束形成方法，首先设定了光速、信号频率、阵元数目、阵元间距、方向角和时间等参数。

然后，根据给定的信号幅度和阵列流形计算波束形成输出。

最后，使用MATLAB的绘图功能将波束形成输出绘制成图形。

Matlab仿真4天线准正交空时编码结合波束形成

clcclear allSNR_dB=0:5:25;len=length(SNR_dB);for i=1:1:lenSNR_dB(i)P(i)=qostbc(SNR_dB(i));endsemilogy(SNR_dB,P,'*')xlabel('SNR/dB');ylabel('BER');grid on%------------------------------------------------------------------------function Pb=qostbc(SNR_dB)%------------------------------------------------------------------------%星座映射采用的是QPSK映射% SNR_dB 接收天线上的接收信噪比,在仿真的过程中把接收信号的功率归一化为1 % Pb是误比特率%-------------------------------------------------------------------------SNR=10^(SNR_dB/10);Cons=sqrt(1/2)*[1+j -1+j -1-j 1-j]; % QPSK星座E=1; %接收端的信号功率No=E/(SNR);%高斯白噪声的功率谱密度Frame=100000;ErrorNum=0;for i=1:1:Frametemp=rand(1);if (temp<0.25)x1=Cons((temp<0.25));elseif (temp<0.5)x1=Cons((temp<0.5)+1);elseif (temp<0.75)x1=Cons((temp<0.75)+2); elsex1=Cons((temp<1)+3); endendendtemp=rand(1);if (temp<0.25)x2=Cons((temp<0.25)); elseif (temp<0.5)x2=Cons((temp<0.5)+1); elseif (temp<0.75)x2=Cons((temp<0.75)+2);elsex2=Cons((temp<1)+3); endendendtemp=rand(1);if (temp<0.25)x3=Cons((temp<0.25)); elseif (temp<0.5)x3=Cons((temp<0.5)+1); elseif (temp<0.75)x3=Cons((temp<0.75)+2); elsex3=Cons((temp<1)+3); endendendtemp=rand(1);if (temp<0.25)x4=Cons((temp<0.25)); elseif (temp<0.5)x4=Cons((temp<0.5)+1); elseif (temp<0.75)x4=Cons((temp<0.75)+2); elsex4=Cons((temp<1)+3); endendendX=[x1 x2 x3 x4;-conj(x2) conj(x1) -conj(x4) conj(x3);-conj(x3) -conj(x4) conj(x1) conj(x2);x4 -x3 -x2 x1]; % 进行空时编码h1=sqrt(1/2)*(randn(1)+j*randn(1)); % 产生信道h2=sqrt(1/2)*(randn(1)+j*randn(1));h3=sqrt(1/2)*(randn(1)+j*randn(1));h4=sqrt(1/2)*(randn(1)+j*randn(1));H=[h1;h2;h3;h4];N=sqrt(No/2)*(randn(4,1)+j*randn(4,1));%接收端的噪声bf1=sqrt(1/16)*[1;exp(j*pi*sin(0.1*pi));exp(j*2*pi*sin(0.1*pi));exp (j*3*pi*sin(0.1*pi));exp(j*4*pi*sin(0.1*pi));exp(j*5*pi*sin(0.1*pi));exp(j*6*pi*sin(0.1*pi));exp(j*7 *pi*sin(0.1*pi))]'*h1*[1;exp(j*pi*sin(0.1*pi));exp(j*2*pi*sin(0.1*pi));exp(j*3*pi*sin(0.1*pi));ex p(j*4*pi*sin(0.1*pi));exp(j*5*pi*sin(0.1*pi));exp(j*6*pi*sin(0.1*pi));exp(j*7*pi*sin(0.1*pi))];bf2=sqrt(1/16)*[1;exp(j*pi*sin(0.2*pi));exp(j*2*pi*sin(0.2*pi));exp (j*3*pi*sin(0.2*pi));exp(j*4*pi*sin(0.2*pi));exp(j*5*pi*sin(0.2*pi));exp(j*6*pi*sin(0.2*pi));exp(j*7 *pi*sin(0.2*pi))]'*h2*[1;exp(j*pi*sin(0.2*pi));exp(j*2*pi*sin(0.2*pi));exp(j*3*pi*sin(0.2*pi));ex p(j*4*pi*sin(0.2*pi));exp(j*5*pi*sin(0.2*pi));exp(j*6*pi*sin(0.2*pi));exp(j*7*pi*sin(0.2*pi))];bf3=sqrt(1/16)*[1;exp(j*pi*sin(0.8*pi));exp(j*2*pi*sin(0.8*pi));exp (j*3*pi*sin(0.8*pi));exp(j*4*pi*sin(0.8*pi));exp(j*5*pi*sin(0.8*pi));exp(j*6*pi*sin(0.8*pi));exp(j*7 *pi*sin(0.8*pi))]'*h3*[1;。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

Capon波束形成
阵列N=16,信号
0-30
θ︒
=，干扰为
160
θ︒=，
219
θ︒
=，
345
θ︒
=，干扰功率分别为：40dB,35dB,50dB。

Capon波束形成后的方向图和功率谱如下：
为了比较接收数据直接估计噪声协方差矩阵和利用干扰+噪声估计协方差矩阵的Capon 波束形成的差异，进行如下仿真：
可以看出利用干扰+噪声估计协方差矩阵的方向图性能较优于接收数据直接估计噪声协方差矩阵的方向图。

代码：
clc;
clear all ;
close all;
ima=sqrt(-1);
element_num=8; %阵元数
d_lamda=1/2; %阵元间距与波长的关系
theta=-90:0.5:90; %范围
theta0=-30; %来波方向
theta1=60; %干扰方向1
theta2=19; %干扰方向2
theta3=45; %干扰方向3
L=1000; %采样单元数
for i=1:L;
amp0=10*randn(1);%信号的幅度随机产生，保证信号之间是不相关的amp1=100*randn(1);%输入阵列的噪声
amp2=sqrt(10^3.5)*randn(1);%输入阵列的噪声
amp3=sqrt(10^5)*randn(1);%输入阵列的噪声
ampn=3;%噪声
x(:,i)=amp0*exp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]')+...
amp1*exp(ima*2*pi*1/2*sin(theta1*pi/180)*[0:element_num-1]')+...
amp2*exp(ima*2*pi*1/2*sin(theta2*pi/180)*[0:element_num-1]')+...
amp3*exp(ima*2*pi*1/2*sin(theta3*pi/180)*[0:element_num-1]')+...
ampn*(randn(element_num,1)+ima*randn(element_num,1));
end
Rx=1/L*x* x';
R=inv(Rx);
steer=exp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]');
w=R*steer/(steer'*R*steer);%Capon最优权矢量
for j=1:length(theta);
a=exp(ima*2*pi*d_lamda*sin(theta(j)*pi/180)*[0:element_num-1]');
f(j)=w'*a;
p(j)=1/(a'*R*a);
end
F=20*log10(abs(f)/(max(abs(f))));
P=20*log10(abs(p)/(max(abs(p))));%此处是功率的对数形式
figure;
% subplot(121)
plot(theta,F),grid on,hold on
plot(theta0,-80:0,'.')
plot(theta1,-80:0,'.')
plot(theta2,-80:0,'.')
plot(theta3,-80:0,'.')
xlabel('theta/o');
ylabel('F/dB');
title('Capon beamforming方向图') % axis([-90 90 -50 0]);
% subplot(122)
figure;
plot(theta,P),grid on
xlabel('theta/o');
ylabel('功率/dB');
title('Capon beamforming功率谱')。