Coming soon!

The Least Mean Squares (LMS) beamformer is a low-complexity beamformer introduced by Bernard Widrow. It is different from every beamformer we have shown so far in two ways: 1) it requires knowing the SOI, or at least part of it (e.g., a synchronization sequence, pilots, etc) and 2) it is iterative, meaning the weights are honed in on over some number of iterations. It works by minimizing the mean square error between the desired signal (the SOI) and the output of the beamformer (i.e., the weights applied to the received samples). The traditional implementation of LMS is to treat each received sample as the next step in the iterative process, by applying the current weights to the single sample and calculating the error. That error is then used to fine-tune the weights, and the process repeats. The LMS beamformer can be used in both analog and digital beamforming. The LMS algorithm is given by the following equation:

In the example Python code below, we simulate an 8-element array with a SOI that is comprised of a repeating Gold code transmitted as BPSK. Gold codes are used in 5G and GPS, and have excellent cross-correlation properties, making them great for synchronization signals. In the simulation we also include two tone interferers, at 60 and -50 degrees. Note that this simulation does not include any time or frequency shift, if it did then we would have to synchronize to the SOI as part of the LMS process (i.e., joint beamforming-synchronization). In the following animation we sweep the SOI angle of arrival and plot the beam pattern LMS created for us after 10k samples. Note how LMS keeps the gain towards the SOI at exactly 0 dB (unless there is an interferer on top), while placing nulls at the interferers.

..image::../_images/doa_lms_animation.gif

..code-block::python

Try changing:code:`theta_soi`, the amount of noise (i.e., the:code:`0.5*n`), and the step size:code:`mu` to see how the LMS algorithm performs.

Circular Arrays

importmatplotlib.pyplotasplt

frommatplotlib.animationimportFuncAnimation

t=np.arange(N)/sample_rate# time vector

theta1=30/180*np.pi# convert to radians

theta2=60/180*np.pi

theta3=-60/180*np.pi

s1=np.exp(-2j*np.pi*d*np.arange(Nr)*np.sin(theta1)).reshape(-1,1)# 8x1

s2=np.exp(-2j*np.pi*d*np.arange(Nr)*np.sin(theta2)).reshape(-1,1)

s3=np.exp(-2j*np.pi*d*np.arange(Nr)*np.sin(theta3)).reshape(-1,1)

ifFalse:

soi=np.exp(2j*np.pi*0.01e6*t).reshape(1,-1)# 1xN

else:

gold_code=np.array([-1,-1,1,-1,-1,1,1,1,1,1,1,1,1,-1,1,1,-1,1,-1,1,-1,1,-1,-1,1,1,1,1,1,-1,1],dtype=complex)

soi=np.repeat(gold_code,soi_samples_per_symbol)# Gold code is 31 bits, so 31*8 = 248 samples

num_sequence_repeats=int(N/soi.shape[0])+1# number of times to repeat the sequence

soi=np.tile(soi,num_sequence_repeats)# repeat the sequence

soi=soi[:N]# trim to N samples

soi=soi.reshape(1,-1)# 1xN

tone2=np.exp(2j*np.pi*0.02e6*t).reshape(1,-1)

tone3=np.exp(2j*np.pi*0.03e6*t).reshape(1,-1)

fig,ax=plt.subplots(subplot_kw={'projection':'polar'})

theta_sois_deg=np.linspace(-90,90,200)

fortheta_soi_degintheta_sois_deg:

t=np.arange(N)/sample_rate# time vector

print(np.var(soi))

print(np.var(tone2))

print(np.var(tone3))

theta_soi=theta_soi_deg/180*np.pi# convert to radians

theta2=60/180*np.pi

theta3=-50/180*np.pi

s1=np.exp(-2j*np.pi*d*np.arange(Nr)*np.sin(theta_soi)).reshape(-1,1)# 8x1

s2=np.exp(-2j*np.pi*d*np.arange(Nr)*np.sin(theta2)).reshape(-1,1)

s3=np.exp(-2j*np.pi*d*np.arange(Nr)*np.sin(theta3)).reshape(-1,1)

r=s1 @soi+s2 @tone2+0.1*s3 @tone3

n=np.random.randn(Nr,N)+1j*np.random.randn(Nr,N)

ifFalse:

soi=np.exp(2j*np.pi*0.01e6*t).reshape(1,-1)# 1xN

else:

gold_code=np.array([-1,1,1,-1,1,1,1,1,-1,-1,-1,1,1,-1,-1,-1,-1,-1,1,1,1,-1,-1,1,1,1,-1,1,1,1,1,1,1,-1,-1,-1,1,1,1,-1,-1,1,1,-1,-1,1,-1,1,-1,-1,1,-1,-1,-1,-1,-1,-1,1,1,-1,1,-1,1,-1,1,1,-1,-1,-1,-1,1,1,1,-1,1,-1,1,1,1,1,1,-1,-1,-1,-1,1,1,1,-1,1,-1,-1,-1,1,1,1,1,-1,1,1,1,-1,1,-1,-1,-1,-1,1,-1,1,1,-1,-1,-1,-1,1,-1,1,1,-1,-1,-1,-1,-1,-1,1,1])

soi=np.repeat(gold_code,soi_samples_per_symbol)# Gold code is 31 bits, so 31*8 = 248 samples

num_sequence_repeats=int(N/soi.shape[0])+1# number of times to repeat the sequence

soi=np.tile(soi,num_sequence_repeats)# repeat the sequence to fill simulated time

soi=soi[:N]# trim to N samples

soi=soi.reshape(1,-1)# 1xN

defw_mvdr(theta,r):

s=np.exp(-2j*np.pi*d*np.arange(r.shape[0])*np.sin(theta))# steering vector in the desired direction theta

s=s.reshape(-1,1)# make into a column vector (size 3x1)

R=np.cov(r)# Calc covariance matrix. gives a Nr x Nr covariance matrix of the samples

Rinv=np.linalg.pinv(R)# 3x3. pseudo-inverse tends to work better than a true inverse

w= (Rinv @s)/(s.conj().T @Rinv @s)# MVDR/Capon equation! numerator is 3x3 * 3x1, denominator is 1x3 * 3x3 * 3x1, resulting in a 3x1 weights vector

tone2=np.exp(2j*np.pi*0.02e6*t).reshape(1,-1)

tone3=np.exp(2j*np.pi*0.03e6*t).reshape(1,-1)

defpower_mvdr(theta,r):

s=np.exp(-2j*np.pi*d*np.arange(r.shape[0])*np.sin(theta))# steering vector in the desired direction theta_i

s=s.reshape(-1,1)# make into a column vector (size 3x1)

R=np.cov(r)

Rinv=np.linalg.pinv(R)# 3x3. pseudo-inverse tends to work better than a true inverse

return1/(s.conj().T @Rinv @s).squeeze()

r=s1 @soi+s2 @tone2+s3 @tone3

n=np.random.randn(Nr,N)+1j*np.random.randn(Nr,N)

theta_scan=np.linspace(-1*np.pi,np.pi,1000)# 1000 different thetas between -180 and +180 degrees

results= []

fortheta_iintheta_scan:

results.append(10*np.log10(power_mvdr(theta_i,r)))# compare to using equation for MVDR power, should match, SHOW MATH OF WHY THIS HAPPENS!

results-=np.max(results)# normalize

w_lms=np.zeros((Nr,1),dtype=np.complex128)# zeros

ifFalse:

fig,ax=plt.subplots(subplot_kw={'projection':'polar'})

ax.plot(theta_scan,results)# MAKE SURE TO USE RADIAN FOR POLAR

ax.set_theta_zero_location('N')# make 0 degrees point up

ax.set_theta_direction(-1)# increase clockwise

ax.set_rlabel_position(30)# Move grid labels away from other labels

ax.set_thetamin(-90)

ax.set_thetamax(90)

plt.show()

error_log= []

foriinrange(N):

r_sample=r[:,i].reshape(-1,1)# 8x1

soi_sample=soi[0,i]# scalar

y=w_lms.conj().T @r_sample# apply the weights (output is a scalar)

y=y.squeeze()# make it a scalar

error_log.append(np.abs(error)**2)

w_lms=w_lms*gamma+mu*np.conj(error)*r_sample# weights are still 8x1

w_lms/=np.linalg.norm(w_lms)# normalize

w_soi=w_mvdr(theta1,r)

w_soi=w_soi.reshape(-1)# make into a row vector

print(w_soi)

w_soi=np.conj(w_soi)

w_padded=np.concatenate((w_soi,np.zeros(N_fft-Nr)))# zero pad to N_fft elements to get more resolution in the FFT

w_fft_dB=10*np.log10(np.abs(np.fft.fftshift(np.fft.fft(w_padded)))**2)# magnitude of fft in dB

w_fft_dB-=np.max(w_fft_dB)# normalize to 0 dB at peak

theta_bins=np.arcsin(np.linspace(-1,1,N_fft))# in radians

fig,ax=plt.subplots(subplot_kw={'projection':'polar'})

ax.plot(theta_bins,w_fft_dB)# MAKE SURE TO USE RADIAN FOR POLAR

ax.set_theta_zero_location('N')# make 0 degrees point up

ax.set_theta_direction(-1)# increase clockwise

ax.set_rlabel_position(55)# Move grid labels away from other labels

ax.set_thetamin(-90)# only show top half

ax.set_thetamax(90)

ax.set_ylim([-30,1])# because there's no noise, only go down 30 dB

plt.show()

w_lms=w_lms.reshape(-1)# make into a row vector

w_lms=w_lms.conj()

w_padded=np.concatenate((w_lms,np.zeros(N_fft-Nr)))# zero pad to N_fft elements to get more resolution in the FFT

w_fft_dB=10*np.log10(np.abs(np.fft.fftshift(np.fft.fft(w_padded)))**2)# magnitude of fft in dB

theta_bins=np.arcsin(np.linspace(-1,1,N_fft))# in radians

ax.plot(theta_bins,w_fft_dB)

ax.plot([theta_soi,theta_soi], [-30,10],'g--')

ax.plot([theta2,theta2], [-30,10],'r--')

ax.plot([theta3,theta3], [-30,10],'r--')

ax.set_theta_zero_location('N')# type: ignore # make 0 degrees point up

ax.set_theta_direction(-1)# type: ignore # increase clockwise

ax.set_rlabel_position(55)# type: ignore # Move grid labels away from other labels

ax.set_thetamin(-90)# type: ignore # only show top half

ax.set_thetamax(90)# type: ignore

ax.set_ylim((-30,10))# because there's no noise, only go down 30 dB

plt.draw()

plt.pause(0.001)# pause to update the plot

plt.cla()