MVDR 的原始公式算法需要对一段内的数据进行求Rx和求逆,原始代码实现的参考如下
%单目标MVDR和CBF波束形成,通过调整SNR改变信噪比,根据实验所需信号来选择性注释
clear;clc;fclose('all');
%声源、水平阵、环境参数设置
%环境参数
c = 1500; %声速1500m/s
SNR = 20; %信噪比为20dB
%声源(假定在很远处的点源,近场球面波,到达很远处的水平阵的时候可近似为平面波)
SL = 0; %信号能量140dB
r = 7000; %与水平阵首阵元相距7000m
f = 50; %信号频率50Hz
lambda=c/f; %波长
angle=10; %入射角
anglenn=70; %入射角
%水平阵
M = 16; %水平阵阵元数9
d = c/(2*f); %阵元间距15m
angles=-90:1:90; %检测角度范围
angles=10; %检测角度范围
Nsnapshot=15000; %快拍数(结合多次发射信号)
%计算各阵元接收信号x=s+n,即信号加噪声
%阵列响应向量
v=1\exp(-1j*2*pi*(d*sin(angle*pi/180)/lambda)*(-(M-1)/2:(M-1)/2)');
vnn=1\exp(-1j*2*pi*(d*sin(anglenn*pi/180)/lambda)*(-(M-1)/2:(M-1)/2)');
%发射信号
ss=zeros(1,Nsnapshot);
ss(1,1:200:Nsnapshot)=1;
s=sqrt(10^(SL/10))*exp(1j*2*pi*f*(1:Nsnapshot));
s=s.*ss;
ssnn=zeros(1,Nsnapshot);
ssnn(1,1:20:Nsnapshot)=1;
snn=sqrt(10^(SL/10))*exp(1j*2*pi*f*(1:Nsnapshot));
snn=snn.*ssnn;
%高斯白噪声
n=sqrt(10^((SL-SNR)/10))*(randn(M,Nsnapshot)+1j*randn(M,Nsnapshot))/sqrt(2);
%接收信号
rx=v*s+vnn*snn+n;
%计算响应向量和波束形成响应
% Rx=(M*v*(s*s')*v'+10^((SL-SNR)/10)*eye(M))/Nsnapshot;%理想信号tr
% A=s.';
A=rx';
Rx=(A'*A)/Nsnapshot;%加入高斯白噪声的真实信号
%驾驶向量
c=1\exp(-1j*2*pi*(-(M-1)/2:(M-1)/2)'*(d*sin(angles*pi/180)/lambda));
%直接求解闭式解
Cmv1=(Rx\c)/(diag(diag((c'/Rx)*c)));
%
%% 使用原始公式进行计算
w=(inv(Rx)*c)/(c'*inv(Rx)*c);
y=w.'*A.';
figure;
plot(abs(y));
figure;
plot(abs(s));
如果使用FPGA进行Verilog实现需要消耗大量资源,因此我们需要对这个原始算法进行修正,我们参考了论文《VLSI_systolic_arrays_for_adaptive_nulling_radar》的实现方式,引入了QR分解和迭代遗忘因子和高斯消元,使得整个计算“相对而言“提高计算速度和节省fpga资源,下面是改进后的思路和方案的matlab代码,
%单目标MVDR和CBF波束形成,通过调整SNR改变信噪比,根据实验所需信号来选择性注释
clear;clc;fclose('all');
%声源、水平阵、环境参数设置
%环境参数
c = 1500; %声速1500m/s
SNR = 20; %信噪比为20dB
%声源(假定在很远处的点源,近场球面波,到达很远处的水平阵的时候可近似为平面波)
SL = 0; %信号能量140dB
r = 7000; %与水平阵首阵元相距7000m
f = 50; %信号频率50Hz
lambda=c/f; %波长
angle=10; %入射角
anglenn=70; %入射角
%水平阵
M = 16; %水平阵阵元数9
d = c/(2*f); %阵元间距15m
angles=-90:1:90; %检测角度范围
angles=10; %检测角度范围
Nsnapshot=15000; %快拍数(结合多次发射信号)
%计算各阵元接收信号x=s+n,即信号加噪声
%阵列响应向量
v=1\exp(-1j*2*pi*(d*sin(angle*pi/180)/lambda)*(-(M-1)/2:(M-1)/2)');
vnn=1\exp(-1j*2*pi*(d*sin(anglenn*pi/180)/lambda)*(-(M-1)/2:(M-1)/2)');
%发射信号
ss=zeros(1,Nsnapshot);
ss(1,1:200:Nsnapshot)=1;
s=sqrt(10^(SL/10))*exp(1j*2*pi*f*(1:Nsnapshot));
s=s.*ss;
ssnn=zeros(1,Nsnapshot);
ssnn(1,1:20:Nsnapshot)=1;
snn=sqrt(10^(SL/10))*exp(1j*2*pi*f*(1:Nsnapshot));
snn=snn.*ssnn;
%高斯白噪声
n=sqrt(10^((SL-SNR)/10))*(randn(M,Nsnapshot)+1j*randn(M,Nsnapshot))/sqrt(2);
%接收信号
rx=v*s+vnn*snn+n;
%计算响应向量和波束形成响应
% Rx=(M*v*(s*s')*v'+10^((SL-SNR)/10)*eye(M))/Nsnapshot;%理想信号tr
% A=s.';
A=rx';
Rx=(A'*A)/Nsnapshot;%加入高斯白噪声的真实信号
%驾驶向量
c=1\exp(-1j*2*pi*(-(M-1)/2:(M-1)/2)'*(d*sin(angles*pi/180)/lambda));
%直接求解闭式解
Cmv1=(Rx\c)/(diag(diag((c'/Rx)*c)));
%
%% 使用原始公式进行计算
w=(inv(Rx)*c)/(c'*inv(Rx)*c);
y=w.'*A.';
figure;
plot(abs(y));
figure;
plot(abs(s));
%% 使用原始公式拆分后的过程计算
x=(inv(Rx)*c);
w=x/(c'*x);
y=w.'*A.';
figure;
plot(abs(y));
figure;
plot(abs(s));
%% 使用qr分解进行简化计算
R=qr(A);
R=R(1:M,:);
Rxqr=(R'*R)/Nsnapshot;
x=(inv(Rxqr)*c);
w=x/(c'*x);
y=w.'*A.';
figure;
plot(abs(y));
figure;
plot(abs(s));
%% 以最后一次迭代结果计算
R= zeros(M);
for i = 1:Nsnapshot
R=qr([R;A(i,:)]);
R=R(1:M,:)*0.99999;
end
Rxx=(R'*R)/Nsnapshot;
x=(inv(Rxx)*c);
w=x/(c'*x);
y=w.'*A.';
figure;
plot(abs(y));
figure;
plot(abs(s));
%% 实时计算
R= zeros(M);
y=[];
for i = 1:Nsnapshot
R=qr([R;A(i,:)]);
R=R(1:M,:)*0.99999;
Rxx=(R'*R)/Nsnapshot;
x=(inv(Rxx)*c);
w=x/(c'*x);
y(i)=w.'*A(i,:).';
end
figure;
plot(abs(y));
figure;
plot(abs(s));
%% 实时计算 使用高斯消元
R= zeros(M);
y=[];
for i = 1:Nsnapshot
R=qr([R;A(i,:)]);
R=R(1:M,:)*0.99999;
% Rxx=(R'*R)/Nsnapshot;
% x=(inv(Rxx)*c);
x=R\(R'\c);
w=x/(c'*x);
y(i)=w.'*A(i,:).';
end
figure;
plot(abs(y));
figure;
plot(abs(s));
具体MVDR Verilog MVDR VHDL MVDR RTL MVDR FPGA代码方案,博主可以进行合作,学生毕设,课设勿扰 , 交流v:183 5645 7685