LMS beamforming Julia Programming
I have been trying to implement a simple LMS adaptive beamforming code. Since I don't have a Matlab license I decided to use Julia since they are quite similar. In order to get a basic code working I implemented the MVRD beamforming example found on Matlabs website (I can't seem o find the link right now). Then I used the link https://teaandtechtime.com/adaptive-beamforming-with-lms/ to get LMS going.
My code at the moment is
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)
The LMS code does not converge (it diverges rather) and I don't know why. I also don't understand the reference signal bit and how it is different from the target steering vector. Perhaps some clarification on the general concepts would be really helpful. I am new to beamforming and my background is in root solvers and such.
Below is the the working Julia code that is rewritten from the Matlab example. It is almost identical to the code above but without the LMS section.
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)
A good practice would be building the simulation iteratively.
One should be aware that the adaptive filter adapts to reference signal while all other data is interfering it. The more correlated the data to the reference the harder is the way to adapt to it.
So the reference signal should be something you can filter out from the signal on the sensors. For example, it can be the non delayed version of it or the signal of one of the sensors gets.
To make things simpler I created an example of real signals which basically gets a different delay according to the sensors.
The simplest way to add a delay to a real signal is by adjusting the phase of its analytic signal and taking the real.
Then we basically have N signals with M samples each which are getting into the adaptive filter which applies an inner product on each row of this M x N matrix.
I created a simple simulation both in MATLAB and Julia.
This is the result from the Julia code:
In the above we see the antenna gain for a target signal coming at 30 [Deg] and interference at [20.0; -35.0; 50.0].
As can be seen the array does adapt and have a gain of ~1 on the reference direction while rejection all other directions.
The full code is available on my StackExchange Signal Processing Q81138 GitHub Repository (Look at the SignalProcessing\Q81138 folder).