波束形成（最大信噪比准则、LCMV、MSC、LMS、RLS）

波束形成的基本概念

# 波束形成准则

分别对上面所阐述的最大信噪比准则，旁瓣对消器，线性最小均方误差，以及自适应LMS和RLS算法进行仿真。
最大信噪比准则：仿真中阵列有10个阵元，阵元间隔为0.5倍波长，信号的入射角为20度，得到的仿真结果的方向图如图1所示。
LCMV:仿真中阵列有8个阵元，阵元间隔为0.5倍波长，信号入射角为0度，干扰入射角为-30和30，得到的仿真结果的方向图如图2所示。
MSC:仿真中阵列有8个阵元，阵元间隔为0.5倍波长，信号入射角为0度，干扰入射角为-40和40，得到的仿真结果的方向图如图3所示。
LMS:仿真中阵列有8个阵元，阵元间隔为0.5倍波长，信号入射角为0度，干扰入射角为-30和30，得到的仿真结果的方向图如图4所示。
RLS:仿真中阵列有8个阵元，阵元间隔为0.5倍波长，信号入射角为0度，干扰入射角为-30和30，得到的仿真结果的方向图如图5所示。

附程序
%lcmv
clear all;
J=8;%传感器的个数
j=sqrt(-1);
r=1/2;%避免空间采样重叠设置传感器间隔与波长的比值为0.5
sita=-90:0.5:90;%设置角度范围
sita1=0;%设置信号入射角
sita2=[-30 30];%设置干扰入射角
N=500;%设置采样点数
s=exp(j2pi*(16000*[0:1:N-1]/200)) ;%构造有用信号（正弦信号）
noise=(randn(J,N)+jrandn(J,N))/sqrt(2);%噪声分量
a1=sqrt(10)exp(-jpi[0:J-1]‘sin(sita1/180pi));
a2=sqrt(10)exp(-jpi*[0:J-1]‘sin(sita2/180pi));
xi=(randn(2,N)+jrandn(2,N))/sqrt(2);%干扰
X=a1s+a2xi+noise;%含噪声的接收到的信号
Rx=XX’/N;
C=a1;
%F=1;
w=inv(Rx)Cinv(C’inv(Rx)C);
for n=1:length(sita)
a=exp(-j2pirsin(sita(n)pi/180)[0:J-1]’);%导向矢量
p(n)=w’a;%天线方向图
end
P=20log10(abs§/max(abs(abs§)));%换算成dB
plot(sita,P);
grid on;
xlabel(‘sita/度’);
ylabel(‘P / dB’);
title(‘lcmv方向图’);

lms程序：
clc;
J=8;%传感器的个数
j=sqrt(-1);
r=1/2;%避免空间采样重叠设置传感器间隔与波长的比值为0.5
sita=-90:0.5:90;%设置角度范围
sita1=0;%设置信号入射角
sita2=[-30 30];%设置干扰入射角
N=500;%设置采样点数
s=exp(j2pi*(16000*[0:1:N-1]/200)) ;%有用信号（正弦信号）
noise=(randn(J,N)+jrandn(J,N))/sqrt(2);%噪声
a1=sqrt(10)exp(-jpi[0:J-1]'sin(sita1/180pi));
a2=sqrt(10)exp(-jpi*[0:J-1]'sin(sita2/180pi));
xi=(randn(2,N)+jrandn(2,N))/sqrt(2);%干扰
X=a1s+a2xi+noise;%含噪声的接收到的信号
xs=a1s;
yd=xs(1,:);%期望信号
mu=0.005;%迭代步长
w=zeros(J,1);%w初始状态
for k =1:N
y(k)=w’X(:,k);
e(k)=yd(k)-y(k);%公式
w=w+muX(:,k)*conj(e(k));%公式
end

for n=1:length(sita)
a=exp(-j2pirsin(sita(n)pi/180)[0:J-1]’);%导向矢量
p(n)=w’a;%天线方向图
end
P=20log10(abs§/max(abs(abs§)));%换算成dB
plot(sita,P);
grid on;
xlabel(‘sita/度’);
ylabel(‘P / dB’);
title(‘lms方向图’);

maxsnr程序：
clc;
J=10;%传感器的个数
j=sqrt(-1);
r=1/2;%避免空间采样重叠设置传感器间隔与波长的比值为0.5
sita=[-90:0.5:90];%设置角度范围
sita1=20;%设置信号入射角
sita2=[-40 40];%设置干扰入射角
N=500;%设置采样点数
s=exp(j2pi*(16000*[0:1:N-1]/200)) ;%有用信号（正弦信号）
noise=(randn(J,N)+jrandn(J,N))/sqrt(2);%噪声
a1=sqrt(10)exp(-jpi[0:J-1]‘sin(sita1/180pi));
% a2=sqrt(10)exp(-jpi*[0:J-1]‘sin(sita2/180pi));
% xi=(randn(2,N)+jrandn(2,N))/sqrt(2);%干扰
% X=a1s+a2xi+noise;%含噪声的接收到的信号
xs=a1;
Rs=xsxs’/N;%信号的自相关矩阵
Rn=noisenoise’/N;%噪声的自相关矩阵
[V,D] = eig(Rs,Rn);%求(Rs,Rn)的广义特征值构成的对角阵D和特征向量构成的矩阵V
[D,ind] = sort(diag(D));%取出D的对角线元素并按照从小到大次序排列，ind为原矩阵D的变动索引
w=V(:,ind(10));%最优权值w为(Rs,Rn)的最大广义特征值所对应的特征向量
for n=1:length(sita)
a=exp(-j2pir*sin(sita(n)pi/180)[0:J-1]’);%导向矢量
p(n)=w’a;%天线方向图
end
P=20log10(abs§/max(abs§));%换算成dB
plot(sita,P);
grid on;
xlabel(‘sita/度’);
ylabel(‘P / dB’);
title(‘最大信噪比方向图’);

msc程序：
clc;
J=8;%传感器的个数
j=sqrt(-1);
r=1/2;%避免空间采样重叠设置传感器间隔与波长的比值为0.5
sita=-90:0.2:90;%设置角度范围
sita1=0;%设置信号入射角
sita2=[-40 40];%设置干扰入射角
N=500;%设置采样点数
s=exp(j2pi*(16000*[0:1:N-1]/200)) ;%有用信号（正弦信号）
noise1=(randn(J,N)+jrandn(J,N))/sqrt(2);%辅助信道噪声
noise2=(randn(1,N)+jrandn(1,N))/sqrt(2);%主信道噪声
a1=sqrt(10)exp(-jpi*[0:J-1]'sin(sita1/180pi));%角sita1上的阵列响应矢量
a2=sqrt(10)exp(-jpi*[0:J-1]'sin(sita2/180pi));%角sita2上的阵列响应矢量
xi=(randn(2,N)+jrandn(2,N))/sqrt(2);%辅助信道干扰
xii=(randn(1,N)+jrandn(1,N))/sqrt(2);%主信道干扰
Xa=a1s+a2xi+noise1;%辅助信道信号
Xmc=(100s+xii+noise2).’;%主信道信号
% xs=a1s;
% yd=xs(1,:);%期望信号
Ra=XaXa’/N;%Xa自相关矩阵
Rma=Xa(conj(Xmc))/N;%Xa和Xmc的互相关矩阵
w=inv(Ra)Rma;
for n=1:length(sita)
a=exp(-j2pirsin(sita(n)pi/180)[0:J-1]’);%导向矢量
p(n)=w’a;%天线方向图
end
P=20log10(abs§/max(abs§));%换算成dB
plot(sita,P);
grid on;
xlabel(‘sita/度’);
ylabel(‘P / dB’);
title(‘MSC方向图’);
RLS程序：
clc;
J=8;%传感器的个数
j=sqrt(-1);
r=1/2;%避免空间采样重叠设置传感器间隔与波长的比值为0.5
sita=-90:0.5:90;%设置角度范围
sita1=0;%设置信号入射角
sita2=[-40 40];%设置干扰入射角
N=500;%设置采样点数
s=exp(j2pi(16000*[0:1:N-1]/200)) ;%有用信号（正弦信号）
noise=(randn(J,N)+jrandn(J,N))/sqrt(2);%噪声
a1=sqrt(10)exp(-jpi[0:J-1]'sin(sita1/180pi));
a2=sqrt(10)exp(-jpi*[0:J-1]'sin(sita2/180pi));
xi=(randn(2,N)+jrandn(2,N))/sqrt(2);%干扰
X=a1s+a2xi+noise;%含噪声的接收到的信号
xs=a1s;
yd=xs(1,:);%期望信号
I=eye(J);
p=0.3I;
lam=1;
w=zeros(J,1);%w初始状态
%RLS算法过程
for k =2:N
L=pX(:,k)/(lam+X(:,k)‘pX(:,k));
e(k)=yd(k)-w’X(:,k);
w=w+Lconj(e(k));
p=(p-LX(:,k)'p)/lam;
end
for n=1:length(sita)
a=exp(-j2pirsin(sita(n)pi/180)[0:J-1]’);%导向矢量
pp(n)=w’a;%天线方向图
end
P=20log10(abs(pp)/max(abs(abs(pp))));%换算成dB
plot(sita,P);
grid on;
xlabel(‘sita/度’);
ylabel(‘P / dB’);
title(‘rls方向图’);