figure-generating-scripts/fm.wav

_templates/patrons.html

_spelling/

Radar-2025-Labs/

* - A random direction

- -14.285 dB

Your results may vary due to the random noise being used to calculate the received samples, which get used to calculate:code:`R`. But the main take-away is that the jammers will be in a null and very low power, the 1 degree off from:code:`dir` will be slightly below 0 dB, but still in the main lobe, and then a random direction is going to be lower than 0 dB but higher than the jammers, and very different every run of the simulation.

Your results may vary due to the random noise being used to calculate the received samples, which get used to calculate:code:`R`. But the main take-away is that the jammers will be in a null and very low power, the 1 degree off from:code:`dir` will be slightly below 0 dB, but still in the main lobe, and then a random direction is going to be lower than 0 dB but higher than the jammers, and very different every run of the simulation. Note that with MVDR you get a gain of 0 dB for the main lobe, but if you were to use the conventional beamformer, you would get:math:`10\log_{10}(Nr)`, so about 12 dB for our 16-element array, showing one of the trade-offs of using MVDR.

Conclusion and References

[1] Mailloux, Robert J. Phased Array Antenna Handbook. Second edition, Artech House, 2005

[2] Van Trees, Harry L. Optimum Array Processing: Part IV of Detection, Estimation, and Modulation Theory. Wiley, 2002.

.. |br|raw:: html

<br>

ax1=fig.add_subplot(gs[0:5,3:5],projection='polar')

polar_plot,=ax1.plot(theta_bins,np.zeros(N_fft))

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

ax1.set_theta_direction(-1)# increase clockwise

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

ax1.set_thetamax(90)

ax1.set_ylim([-30,10])

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

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

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

ax1.set_thetamax(90)# type: ignore

ax1.set_ylim((-30,10))

ax1.set_yticks(np.arange(-30,11,10))# Only label every 10 dB

ax2=fig.add_subplot(gs[5:9,3:5])

rect_plot,=ax2.plot(theta_bins*180/np.pi,np.zeros(N_fft))

plt.axis([-90,90,-40,10])

plt.axis((-90,90,-40,10))

plt.xlabel('Theta [Degrees]')

plt.grid()

importmatplotlib.pyplotasplt

theta=np.deg2rad(20)

s=np.exp(-2j*np.pi*d*np.arange(Nr)*np.sin(theta)).reshape(-1,1)# steering vector

t=np.arange(N)/sample_rate

tx=np.exp(2j*np.pi*f_tone*t).reshape(1,-1)# 1 x 10000

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

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

r=s @tx+0.05*n# Nr x 10000

theta_scan=np.linspace(-1*np.pi,np.pi,1000)

results_conventional=np.zeros(len(theta_scan))

results_mvdr=np.zeros(len(theta_scan))

fortheta_iinrange(len(theta_scan)):

s=np.exp(-2j*np.pi*d*np.arange(Nr)*np.sin(theta_scan[theta_i]))

r_weighted_conv=np.conj(w_conventional) @r# apply our weights corresponding to the direction theta_i

results_conventional[theta_i]=10*np.log10(np.mean(np.var(r_weighted_conv)))# energy detector

R= (r @r.conj().T)/r.shape[1]# 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

results_mvdr[theta_i]=10*np.log10(1/(s.conj().T @Rinv @s).squeeze())

phase_diff=np.mean(np.angle(r[0, :]*np.conj(r[1, :])))

theta_est=np.arcsin(phase_diff/ (2*np.pi*d))

print("Estimated angle of arrival: ",theta_est*180/np.pi)# in degrees

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

ax.plot(theta_scan,results_conventional)

ax.plot(theta_scan,results_mvdr)

ax.plot([theta_est]*2, [-20,10],'c')

ax.plot([theta]*2, [-20,10],'g--')

ax.legend(['Conventional','MVDR','np.angle(x x.conj)','Correct AoA'])

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

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

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

ax.set_thetamax(90)# type: ignore

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

plt.show()

exit()

dir=np.asmatrix([np.cos(theta),# x component, goes back to the geometry shown in the beginning of DOA chapter

0]# z component

).T

).T

print("dir:\n",dir)# Remember that it's a unit vector representing a direction, it's not in meters

exit()

theta=np.deg2rad(60)# azimith angle

phi=np.deg2rad(30)# elevation angle

theta=np.deg2rad(0)# azimith angle

phi=np.deg2rad(0)# elevation angle

defget_unit_vector(theta,phi):

returnnp.asmatrix([np.sin(theta)*np.sin(phi),# x component

np.cos(theta)*np.sin(phi),# y component

).T

dir=get_unit_vector(theta,phi)

print("dir:\n",dir)# Remember that it's a unit vector representing a direction, it's not in meters

ifFalse:

az_range_rad=np.radians(np.arange(AZ_MIN,AZ_MAX,AZ_RES))[None,:]#1x360

el_range_rad=np.radians(np.arange(EL_MIN,EL_MAX,EL_RES))[:,None]#180x1

k=2*np.pi*np.asarray([fc]*Nr)/3e8#3x0,

kx=k[0]*np.sin(el_range_rad)*np.cos(az_range_rad)#180*360

ky=k[0]*np.sin(el_range_rad)*np.sin(az_range_rad)#180*360

kz=k[0]*np.cos(el_range_rad)#180x1

print('kx shape:',kx.shape)

print('ky shape:',ky.shape)

print('kz shape:',kz.shape)

scan_kx=np.exp(-1j* (kx[:,:,None]*np.arange(4)*d))[:,:,None,:]# 180x360x1x10

scan_kz=np.exp(-1j* (kz[:,:,None]*np.arange(4)*d))[:, :,:,None]# 180x360x10x1

print('scan_kx shape:',scan_kx.shape)

print('scan_kz shape:',scan_kz.shape)

print('scan_overall shape:',scan_overall.shape)

scan_overall=scan_overall.reshape((scan_overall.shape[0],scan_overall.shape[1],-1))#180x360x100

print('scan_overall shape:',scan_overall.shape)# (180, 360, 16)

weight_conj_t=w.conj().T

weight_conj_t=np.expand_dims(weight_conj_t,axis=0)# 1x1x16

print('weight vector shape:',weight_conj_t.shape)# (1, 16)

results=np.sum(weight_conj_t*scan_overall,axis=2)#180x360

results=np.abs(results)

plt.figure(figsize=(8,6))

im=plt.imshow(results,

extent=(AZ_MIN,AZ_MAX,EL_MIN,EL_MAX),

cmap='viridis',

origin='lower')

plt.colorbar(im)

plt.xlabel('Azimuth Angle (degrees)')

plt.ylabel('Elevation Angle (degrees)')

plt.grid()

plt.show()

ifTrue:

theta_scan=np.linspace(-np.pi,np.pi,resolution)# azimuth angles

phi_scan=np.linspace(0,np.pi,resolution)# elevation angles

results=np.zeros((resolution,resolution))# 2D array to store results

fori,theta_iinenumerate(theta_scan):

forj,phi_iinenumerate(phi_scan):

a=steering_vector(pos,get_unit_vector(theta_i,phi_i))# array factor

results[i,j]=np.abs(w.conj().T @a)[0,0]# power in signal, in dB

plt.imshow(results,extent=(-180,180,0,90),origin='lower',aspect='auto',cmap='viridis')

plt.colorbar(label='Power [dB]')

plt.xlabel('Azimuth angle [degrees]')

plt.ylabel('Elevation angle [degrees]')

plt.show()

ifFalse:

results[i,j]=10*np.log10(np.abs(resp)[0,0])# power in signal, in dB

results[results<-10]=-10# crop the z axis to -10 dB

print(np.max(results))# 12.04 dB because we have 16 elements 10*log10(16) = 12.04 dB

fig,ax=plt.subplots(subplot_kw={"projection":"3d","computed_zorder":False})

surf=ax.plot_surface(np.sin(theta_scan[:,None])*np.sin(phi_scan[None,:]),# x # type: ignore

np.cos(theta_scan[:,None])*np.sin(phi_scan[None,:]),# y

theta_scan=np.linspace(0,2*np.pi,resolution)# azimuth angles

phi_scan=np.linspace(0,np.pi,resolution)# elevation angles

results=np.zeros((resolution,resolution))# 2D array to store results

fori,theta_iinenumerate(theta_scan):

forj,phi_iinenumerate(phi_scan):

dir_i=get_unit_vector(theta_i,phi_i)

a=steering_vector(pos,dir_i)# array factor

resp=w.conj().T @a# scalar

results[i,j]=10*np.log10(np.abs(resp)[0,0])# power in signal, in dB

print(np.max(results))# 1.3 dB

max_idx=np.unravel_index(np.argmax(results,axis=None),results.shape)

print(theta_scan[max_idx[0]]*180/np.pi)# theta doesnt really mean anything here because phi is close to 0

print(phi_scan[max_idx[1]]*180/np.pi)# phi is sometimes +10 and sometimes -10 deg

importmatplotlib.pyplotasplt

frommatplotlib.animationimportFuncAnimation

B_theta=np.ones(Nr)# desired beam pattern

phi= (np.arange(Nr)- (Nr-1)/2)*2*np.pi/Nr# eq 3.92

B=np.conj(B_theta)*np.exp(-1j*phi* (Nr-1)/2)# eq 3.94

b=np.fft.ifft(B)

n=np.arange(Nr)

w=b*np.exp(-1j*n*np.pi* (Nr-1)/Nr)# eq 3.101

w=w/np.sum(w)# normalize weights

w=np.conj(w)# or else our answer will be negative/inverted

w=w.squeeze()

w_padded=np.concatenate((w,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))# Map the FFT bins to angles in radians

fig,ax=plt.subplots()

ax.plot(theta_bins*180/np.pi,w_fft_dB)# MAKE SURE TO USE RADIAN FOR POLAR

ax.set_xlabel("Theta [Degrees]")

ax.set_ylabel("Beam Pattern [dB]")

ax.set_ylim((-30,0))

ax.grid()

plt.show()

importmatplotlib.pyplotasplt

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

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

filenames= []

fortheta_soi_degintheta_sois_deg:

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

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

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

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_lms+=mu*np.conj(error)*r_sample# weights are still 8x1

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

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)

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

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

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()

filename=tmp+'\\lms_animation'+str(theta_soi_deg)+'.png'

print(theta_soi_deg)

fig.savefig(filename,bbox_inches='tight')

filenames.append(filename)

plt.close(fig)

images= []

forfilenameinfilenames:

images.append(imageio.imread(filename))

imageio.mimsave(tmp+'\\lms_animation.gif',images,fps=10,loop=0)# loop=0 means loop forever