LMS 波束成形 Julia 编程
LMS beamforming Julia Programming
我一直在尝试实现一个简单的 LMS 自适应波束形成代码。由于我没有 Matlab 许可证，我决定使用 Julia，因为它们非常相似。为了让基本代码正常工作，我实现了在 Matlabs 网站上找到的 MVRD 波束形成示例（我现在似乎找不到 link）。然后我使用 link https://teaandtechtime.com/adaptive-beamforming-with-lms/ 来启动 LMS。
我现在的密码是
using Plots
using LinearAlgebra
# Source: https://teaandtechtime.com/adaptive-beamforming-with-lms/
M = 20; # Number of Array Elements.
N = 200; # Number of Signal Samples.
n = 1:N; # Time Sample Index Vector.
c = 3*10^8; # Speed of light
f = 2.4*10^9; # Frequency [Hz]
lambda = c/f; # Incoming Signal Wavelength in [m].
d = lambda/2; # Interelement Distance in [m].
SNR = 20; # Target SNR in dBs.
phi_s = 0; # Target azimuth angle in degrees.
phi_i1 = 20; # Interference angle in degrees.
phi_i2 = -30; # Interference angle in degrees.
phi_i3 = 50; # Interference angle in degrees.
INR1 = 35; # Interference #1 INR in dBs.
INR2 = 70; # Interference #2 INR in dBs.
INR3 = 50; # Interference #3 INR in dBs.
u_s = (d/lambda)*sin(phi_s*pi/180); # Normalized Spatial Frequency of the Target signal.
u_int1 = (d/lambda)*sin(phi_i1*pi/180); # Normalized Spatial Frequency of the Interferer #1.
u_int2 = (d/lambda)*sin(phi_i2*pi/180); # Normalized Spatial Frequency of the Interferer #2.
u_int3 = (d/lambda)*sin(phi_i3*pi/180); # Normalized Spatial Frequency of the Interferer #3.
tau_s = (d/c)*sin(phi_s*pi/180); # Time delay of the Target signal.
tau1 = (d/c)*sin(phi_i1*pi/180); # Time delay of the Interferer #1.
tau2 = (d/c)*sin(phi_i2*pi/180); # Time delay of the Interferer #2.
tau3 = (d/c)*sin(phi_i3*pi/180); # Time delay of the Interferer #3.
# Target Signal definition.
s = zeros(ComplexF64,M,N)
v_s = exp.(-1im*2*pi*u_s*collect(0:M-1))/sqrt(M); # Target Steering Vector.
for k=1:N
s[:,k] = 10^(SNR/20)*v_s; # Amplitude of Target Signal Generation.
end
# The uncorrelated unit power thermal noise samples with a Gaussian
# distribution are generated by:
w = (randn(M,N)+1im*randn(M,N))/sqrt(2)
# The interference [1iammer] vectors are generated by:
v_i1 = exp.(-1im*2*pi*u_int1*collect(0:M-1))/sqrt(M)
i_x1 = 10^(INR1/20)*v_i1*(randn(1,N)+1im*randn(1,N))/sqrt(2)
v_i2 = exp.(-1im*2*pi*u_int2*collect(0:M-1))/sqrt(M)
i_x2 = 10^(INR2/20)*v_i2*(randn(1,N)+1im*randn(1,N))/sqrt(2)
v_i3 = exp.(-1im*2*pi*u_int3*collect(0:M-1))/sqrt(M)
i_x3 = 10^(INR3/20)*v_i3*(randn(1,N)+1im*randn(1,N))/sqrt(2)
# The three signals are added to produce the overall array signal.
x = s + i_x1 + i_x2 + i_x3 + w
# Run LMS algorithm
mu = 0.001; # LMS step size
a = ones(ComplexF64,M); # Complex weights
S = zeros(ComplexF64,M); # Complex weights
ts = 1/(N*f); # sample time
for k = 1:N
d = cos(2*pi*f*k*ts); # Reference Signal
S = a.*x[:,k];
y = sum(S);
global e = conj(d) - y;
println(e)
global a += mu*x[:,k]*e; # next weight calculation
end
println(a)
# Array Response
Nsamples1 = 3*10^4
angle1 = -90:180/Nsamples1:90-180/Nsamples1
LMS_Beam_Pat = zeros(ComplexF64,Nsamples1)
for k = 1:Nsamples1
u = (d/lambda)*sin(angle1[k]*pi/180)
v = exp.(-1im*2*pi*u*collect(0:M-1))/sqrt(M); # Azimuth Scanning Steering Vector.
LMS_Beam_Pat[k] = a'*v;
end
# Plot the corresponding Beampatterns.
display(plot(angle1,10*log10.(abs.(LMS_Beam_Pat).^2),xlims=(-90,90),ylims=(-100,0)))
sleep(10)
# PolardB plot
display(plot(angle1*pi/180,10*log10.(abs.(LMS_Beam_Pat).^2), proj=:polar, lims=(-80,0)))
sleep(10)
LMS 代码不收敛（而是发散），我不知道为什么。我也不明白参考信号位以及它与目标导向矢量有何不同。也许对一般概念的一些澄清真的很有帮助。我是波束成形的新手，我的背景是根求解器等。
下面是从 Matlab 示例重写的有效 Julia 代码。它与上面的代码几乎相同，但没有 LMS 部分。
using Plots
using LinearAlgebra
M = 20; # Number of Array Elements.
N = 200; # Number of Signal Samples.
n = 1:N; # Time Sample Index Vector.
c = 3*10^8; # Speed of light
f = 2.4*10^9; # Frequency [Hz]
lambda = c/f; # Incoming Signal Wavelength in [m].
d = lambda/2; # Interelement Distance in [m].
SNR = 20; # Target SNR in dBs.
phi_s = 0; # Target azimuth angle in degrees.
phi_i1 = 20; # Interference angle in degrees.
phi_i2 = -30; # Interference angle in degrees.
phi_i3 = 50; # Interference angle in degrees.
INR1 = 35; # Interference #1 INR in dBs.
INR2 = 70; # Interference #2 INR in dBs.
INR3 = 50; # Interference #3 INR in dBs.
u_s = (d/lambda)*sin(phi_s*pi/180); # Normalized Spatial Frequency of the Target signal.
u_int1 = (d/lambda)*sin(phi_i1*pi/180); # Normalized Spatial Frequency of the Interferer #1.
u_int2 = (d/lambda)*sin(phi_i2*pi/180); # Normalized Spatial Frequency of the Interferer #2.
u_int3 = (d/lambda)*sin(phi_i3*pi/180); # Normalized Spatial Frequency of the Interferer #3.
tau_s = (d/c)*sin(phi_s*pi/180); # Time delay of the Target signal.
tau1 = (d/c)*sin(phi_i1*pi/180); # Time delay of the Interferer #1.
tau2 = (d/c)*sin(phi_i2*pi/180); # Time delay of the Interferer #2.
tau3 = (d/c)*sin(phi_i3*pi/180); # Time delay of the Interferer #3.
# Target Signal definition.
s = zeros(ComplexF64,M,N)
v_s = exp.(-1im*2*pi*u_s*collect(0:M-1))/sqrt(M); # Target Steering Vector.
for k=1:N
s[:,k] = 10^(SNR/20)*v_s; # Amplitude of Target Signal Generation.
end
# The uncorrelated unit power thermal noise samples with a Gaussian
# distribution are generated by:
w = (randn(M,N)+1im*randn(M,N))/sqrt(2)
# The interference [1iammer] vectors are generated by:
v_i1 = exp.(-1im*2*pi*u_int1*collect(0:M-1))/sqrt(M)
i_x1 = 10^(INR1/20)*v_i1*(randn(1,N)+1im*randn(1,N))/sqrt(2)
v_i2 = exp.(-1im*2*pi*u_int2*collect(0:M-1))/sqrt(M)
i_x2 = 10^(INR2/20)*v_i2*(randn(1,N)+1im*randn(1,N))/sqrt(2)
v_i3 = exp.(-1im*2*pi*u_int3*collect(0:M-1))/sqrt(M)
i_x3 = 10^(INR3/20)*v_i3*(randn(1,N)+1im*randn(1,N))/sqrt(2)
#The three signals are added to produce the overall array signal.
x = s + i_x1 + i_x2 + i_x3 + w
iplusn = i_x1 + i_x2 + i_x3 + w
# Calculation of the i+n autocorrelation matrix.
R_ipn = 10^(INR1/10)*(v_i1*v_i1') + 10^(INR2/10)*(v_i2*v_i2') + 10^(INR3/10)*(v_i3*v_i3') + I
InvR = inv(R_ipn)
# Calculate the Beam Patterns.
# MVDR Optimum Beamformer computed for a phi_s = 0 deg.
c_opt = InvR*v_s/(v_s'*InvR*v_s);
# Spatial Matched Filter | Steering Vector Beamformer Eq. (11.2.16).
c_mf = v_s
Nsamples1 = 3*10^4
angle1 = -90:180/Nsamples1:90-180/Nsamples1
Opt_Beam_Pat = zeros(ComplexF64,Nsamples1)
Conv_Beam_Pat = zeros(ComplexF64,Nsamples1)
for k = 1:Nsamples1
u = (d/lambda)*sin(angle1[k]*pi/180)
v = exp.(-1im*2*pi*u*collect(0:M-1))/sqrt(M); # Azimuth Scanning Steering Vector.
Opt_Beam_Pat[k] = c_opt'*v
Conv_Beam_Pat[k] = c_mf'*v
end
# Plot the corresponding Beampatterns.
plot(angle1,10*log10.(abs.(Conv_Beam_Pat).^2))
display(plot!(angle1,10*log10.(abs.(Opt_Beam_Pat).^2),xlims=(-90,90),ylims=(-100,0)))
sleep(10)
# PolardB plot
display(plot(angle1*pi/180,10*log10.(abs.(Opt_Beam_Pat).^2), proj=:polar, lims=(-80,0)))
sleep(10)
# Calculate the SINR loss factor for the Optimum Beamformer:
Nsamples = 3*10^4; # The resolution must be very fine to reveal the true depth of the notches.
Lsinr_opt = zeros(ComplexF64,Nsamples,1);
Lsinr_mf = zeros(Nsamples,1);
SNR0 = M*10^(SNR/10);
angle = -90:180/Nsamples:90-180/Nsamples;
for k=1:Nsamples
v = exp.(-1im*pi*collect(0:M-1)*sin(angle[k]*pi/180))/sqrt(M); # Azimuth Scanning Steering Vector.
c_mf = v; # This is the spatial matched filter beamformer.
Lsinr_opt[k] = v'*InvR*v;
SINRout_mf = real(M*(10^(SNR/10))*(abs(c_mf'*v)^2)/(c_mf'*R_ipn*c_mf));
Lsinr_mf[k] = SINRout_mf/SNR0;
end
plot(angle,10*log10.(abs.(Lsinr_opt)),xlims=(-90,0));
display(plot!(angle,10*log10.(abs.(Lsinr_mf)),xlims=(-90,90),ylims=(-75,5)));
sleep(10)
一个好的做法是迭代地构建模拟。
人们应该意识到，自适应滤波器适应参考信号，而所有其他数据都在干扰它。数据与参考的相关性越高，适应它的方式就越难。
因此，参考信号应该是您可以从传感器信号中过滤掉的信号。例如，它可以是它的非延迟版本或其中一个传感器的信号。
为了让事情更简单，我创建了一个真实信号的示例，它基本上根据传感器获得不同的延迟。
向真实信号添加延迟的最简单方法是调整其解析信号的相位并取真实信号。
然后我们基本上有 N 个信号，每个信号都有 M 个样本，每个样本都进入自适应滤波器，自适应滤波器在这个 M x N 矩阵的每一行上应用内积。
我在 MATLAB 和 Julia 中创建了一个简单的模拟。
这是 Julia 代码的结果：
在上面我们看到目标信号的天线增益在 30 [Deg] 和干扰在 [20.0; -35.0; 50.0]。
可以看出，阵列确实适应并在参考方向上具有 ~1 的增益，同时拒绝所有其他方向。
我的 StackExchange Signal Processing Q81138 GitHub Repository 上有完整的代码（查看 SignalProcessing\Q81138 文件夹）。
我一直在尝试实现一个简单的 LMS 自适应波束形成代码。由于我没有 Matlab 许可证，我决定使用 Julia，因为它们非常相似。为了让基本代码正常工作，我实现了在 Matlabs 网站上找到的 MVRD 波束形成示例（我现在似乎找不到 link）。然后我使用 link https://teaandtechtime.com/adaptive-beamforming-with-lms/ 来启动 LMS。
我现在的密码是
using Plots
using LinearAlgebra
# Source: https://teaandtechtime.com/adaptive-beamforming-with-lms/
M = 20; # Number of Array Elements.
N = 200; # Number of Signal Samples.
n = 1:N; # Time Sample Index Vector.
c = 3*10^8; # Speed of light
f = 2.4*10^9; # Frequency [Hz]
lambda = c/f; # Incoming Signal Wavelength in [m].
d = lambda/2; # Interelement Distance in [m].
SNR = 20; # Target SNR in dBs.
phi_s = 0; # Target azimuth angle in degrees.
phi_i1 = 20; # Interference angle in degrees.
phi_i2 = -30; # Interference angle in degrees.
phi_i3 = 50; # Interference angle in degrees.
INR1 = 35; # Interference #1 INR in dBs.
INR2 = 70; # Interference #2 INR in dBs.
INR3 = 50; # Interference #3 INR in dBs.
u_s = (d/lambda)*sin(phi_s*pi/180); # Normalized Spatial Frequency of the Target signal.
u_int1 = (d/lambda)*sin(phi_i1*pi/180); # Normalized Spatial Frequency of the Interferer #1.
u_int2 = (d/lambda)*sin(phi_i2*pi/180); # Normalized Spatial Frequency of the Interferer #2.
u_int3 = (d/lambda)*sin(phi_i3*pi/180); # Normalized Spatial Frequency of the Interferer #3.
tau_s = (d/c)*sin(phi_s*pi/180); # Time delay of the Target signal.
tau1 = (d/c)*sin(phi_i1*pi/180); # Time delay of the Interferer #1.
tau2 = (d/c)*sin(phi_i2*pi/180); # Time delay of the Interferer #2.
tau3 = (d/c)*sin(phi_i3*pi/180); # Time delay of the Interferer #3.
# Target Signal definition.
s = zeros(ComplexF64,M,N)
v_s = exp.(-1im*2*pi*u_s*collect(0:M-1))/sqrt(M); # Target Steering Vector.
for k=1:N
s[:,k] = 10^(SNR/20)*v_s; # Amplitude of Target Signal Generation.
end
# The uncorrelated unit power thermal noise samples with a Gaussian
# distribution are generated by:
w = (randn(M,N)+1im*randn(M,N))/sqrt(2)
# The interference [1iammer] vectors are generated by:
v_i1 = exp.(-1im*2*pi*u_int1*collect(0:M-1))/sqrt(M)
i_x1 = 10^(INR1/20)*v_i1*(randn(1,N)+1im*randn(1,N))/sqrt(2)
v_i2 = exp.(-1im*2*pi*u_int2*collect(0:M-1))/sqrt(M)
i_x2 = 10^(INR2/20)*v_i2*(randn(1,N)+1im*randn(1,N))/sqrt(2)
v_i3 = exp.(-1im*2*pi*u_int3*collect(0:M-1))/sqrt(M)
i_x3 = 10^(INR3/20)*v_i3*(randn(1,N)+1im*randn(1,N))/sqrt(2)
# The three signals are added to produce the overall array signal.
x = s + i_x1 + i_x2 + i_x3 + w
# Run LMS algorithm
mu = 0.001; # LMS step size
a = ones(ComplexF64,M); # Complex weights
S = zeros(ComplexF64,M); # Complex weights
ts = 1/(N*f); # sample time
for k = 1:N
d = cos(2*pi*f*k*ts); # Reference Signal
S = a.*x[:,k];
y = sum(S);
global e = conj(d) - y;
println(e)
global a += mu*x[:,k]*e; # next weight calculation
end
println(a)
# Array Response
Nsamples1 = 3*10^4
angle1 = -90:180/Nsamples1:90-180/Nsamples1
LMS_Beam_Pat = zeros(ComplexF64,Nsamples1)
for k = 1:Nsamples1
u = (d/lambda)*sin(angle1[k]*pi/180)
v = exp.(-1im*2*pi*u*collect(0:M-1))/sqrt(M); # Azimuth Scanning Steering Vector.
LMS_Beam_Pat[k] = a'*v;
end
# Plot the corresponding Beampatterns.
display(plot(angle1,10*log10.(abs.(LMS_Beam_Pat).^2),xlims=(-90,90),ylims=(-100,0)))
sleep(10)
# PolardB plot
display(plot(angle1*pi/180,10*log10.(abs.(LMS_Beam_Pat).^2), proj=:polar, lims=(-80,0)))
sleep(10)
LMS 代码不收敛（而是发散），我不知道为什么。我也不明白参考信号位以及它与目标导向矢量有何不同。也许对一般概念的一些澄清真的很有帮助。我是波束成形的新手，我的背景是根求解器等。
下面是从 Matlab 示例重写的有效 Julia 代码。它与上面的代码几乎相同，但没有 LMS 部分。
using Plots
using LinearAlgebra
M = 20; # Number of Array Elements.
N = 200; # Number of Signal Samples.
n = 1:N; # Time Sample Index Vector.
c = 3*10^8; # Speed of light
f = 2.4*10^9; # Frequency [Hz]
lambda = c/f; # Incoming Signal Wavelength in [m].
d = lambda/2; # Interelement Distance in [m].
SNR = 20; # Target SNR in dBs.
phi_s = 0; # Target azimuth angle in degrees.
phi_i1 = 20; # Interference angle in degrees.
phi_i2 = -30; # Interference angle in degrees.
phi_i3 = 50; # Interference angle in degrees.
INR1 = 35; # Interference #1 INR in dBs.
INR2 = 70; # Interference #2 INR in dBs.
INR3 = 50; # Interference #3 INR in dBs.
u_s = (d/lambda)*sin(phi_s*pi/180); # Normalized Spatial Frequency of the Target signal.
u_int1 = (d/lambda)*sin(phi_i1*pi/180); # Normalized Spatial Frequency of the Interferer #1.
u_int2 = (d/lambda)*sin(phi_i2*pi/180); # Normalized Spatial Frequency of the Interferer #2.
u_int3 = (d/lambda)*sin(phi_i3*pi/180); # Normalized Spatial Frequency of the Interferer #3.
tau_s = (d/c)*sin(phi_s*pi/180); # Time delay of the Target signal.
tau1 = (d/c)*sin(phi_i1*pi/180); # Time delay of the Interferer #1.
tau2 = (d/c)*sin(phi_i2*pi/180); # Time delay of the Interferer #2.
tau3 = (d/c)*sin(phi_i3*pi/180); # Time delay of the Interferer #3.
# Target Signal definition.
s = zeros(ComplexF64,M,N)
v_s = exp.(-1im*2*pi*u_s*collect(0:M-1))/sqrt(M); # Target Steering Vector.
for k=1:N
s[:,k] = 10^(SNR/20)*v_s; # Amplitude of Target Signal Generation.
end
# The uncorrelated unit power thermal noise samples with a Gaussian
# distribution are generated by:
w = (randn(M,N)+1im*randn(M,N))/sqrt(2)
# The interference [1iammer] vectors are generated by:
v_i1 = exp.(-1im*2*pi*u_int1*collect(0:M-1))/sqrt(M)
i_x1 = 10^(INR1/20)*v_i1*(randn(1,N)+1im*randn(1,N))/sqrt(2)
v_i2 = exp.(-1im*2*pi*u_int2*collect(0:M-1))/sqrt(M)
i_x2 = 10^(INR2/20)*v_i2*(randn(1,N)+1im*randn(1,N))/sqrt(2)
v_i3 = exp.(-1im*2*pi*u_int3*collect(0:M-1))/sqrt(M)
i_x3 = 10^(INR3/20)*v_i3*(randn(1,N)+1im*randn(1,N))/sqrt(2)
#The three signals are added to produce the overall array signal.
x = s + i_x1 + i_x2 + i_x3 + w
iplusn = i_x1 + i_x2 + i_x3 + w
# Calculation of the i+n autocorrelation matrix.
R_ipn = 10^(INR1/10)*(v_i1*v_i1') + 10^(INR2/10)*(v_i2*v_i2') + 10^(INR3/10)*(v_i3*v_i3') + I
InvR = inv(R_ipn)
# Calculate the Beam Patterns.
# MVDR Optimum Beamformer computed for a phi_s = 0 deg.
c_opt = InvR*v_s/(v_s'*InvR*v_s);
# Spatial Matched Filter | Steering Vector Beamformer Eq. (11.2.16).
c_mf = v_s
Nsamples1 = 3*10^4
angle1 = -90:180/Nsamples1:90-180/Nsamples1
Opt_Beam_Pat = zeros(ComplexF64,Nsamples1)
Conv_Beam_Pat = zeros(ComplexF64,Nsamples1)
for k = 1:Nsamples1
u = (d/lambda)*sin(angle1[k]*pi/180)
v = exp.(-1im*2*pi*u*collect(0:M-1))/sqrt(M); # Azimuth Scanning Steering Vector.
Opt_Beam_Pat[k] = c_opt'*v
Conv_Beam_Pat[k] = c_mf'*v
end
# Plot the corresponding Beampatterns.
plot(angle1,10*log10.(abs.(Conv_Beam_Pat).^2))
display(plot!(angle1,10*log10.(abs.(Opt_Beam_Pat).^2),xlims=(-90,90),ylims=(-100,0)))
sleep(10)
# PolardB plot
display(plot(angle1*pi/180,10*log10.(abs.(Opt_Beam_Pat).^2), proj=:polar, lims=(-80,0)))
sleep(10)
# Calculate the SINR loss factor for the Optimum Beamformer:
Nsamples = 3*10^4; # The resolution must be very fine to reveal the true depth of the notches.
Lsinr_opt = zeros(ComplexF64,Nsamples,1);
Lsinr_mf = zeros(Nsamples,1);
SNR0 = M*10^(SNR/10);
angle = -90:180/Nsamples:90-180/Nsamples;
for k=1:Nsamples
v = exp.(-1im*pi*collect(0:M-1)*sin(angle[k]*pi/180))/sqrt(M); # Azimuth Scanning Steering Vector.
c_mf = v; # This is the spatial matched filter beamformer.
Lsinr_opt[k] = v'*InvR*v;
SINRout_mf = real(M*(10^(SNR/10))*(abs(c_mf'*v)^2)/(c_mf'*R_ipn*c_mf));
Lsinr_mf[k] = SINRout_mf/SNR0;
end
plot(angle,10*log10.(abs.(Lsinr_opt)),xlims=(-90,0));
display(plot!(angle,10*log10.(abs.(Lsinr_mf)),xlims=(-90,90),ylims=(-75,5)));
sleep(10)
一个好的做法是迭代地构建模拟。
人们应该意识到，自适应滤波器适应参考信号，而所有其他数据都在干扰它。数据与参考的相关性越高，适应它的方式就越难。
因此，参考信号应该是您可以从传感器信号中过滤掉的信号。例如，它可以是它的非延迟版本或其中一个传感器的信号。
为了让事情更简单，我创建了一个真实信号的示例，它基本上根据传感器获得不同的延迟。
向真实信号添加延迟的最简单方法是调整其解析信号的相位并取真实信号。
然后我们基本上有 N 个信号，每个信号都有 M 个样本，每个样本都进入自适应滤波器，自适应滤波器在这个 M x N 矩阵的每一行上应用内积。
我在 MATLAB 和 Julia 中创建了一个简单的模拟。
这是 Julia 代码的结果：
在上面我们看到目标信号的天线增益在 30 [Deg] 和干扰在 [20.0; -35.0; 50.0]。
可以看出，阵列确实适应并在参考方向上具有 ~1 的增益，同时拒绝所有其他方向。
我的 StackExchange Signal Processing Q81138 GitHub Repository 上有完整的代码（查看 SignalProcessing\Q81138 文件夹）。