波达方向估计(DOA)-Python代码实现Beamformer
模拟一个由三根全向天线组成的阵列，然后使用数组来模拟到达阵列的信号。相邻天线之间距离为1/2波长(也称为“半波长间隔”)。将模拟发射机的信号以一定角度theta到达该阵列。另外在这个接收到的信号中添加噪声。
import numpy as np
import matplotlib.pyplot as plt
sample_rate = 1e6
N = 10000 # number of samples to simulate
# Create a tone to act as the transmitted signal
t = np.arange(N) / sample_rate
f_tone = 0.02e6
tx = np.exp(2j * np.pi * f_tone * t)
# Simulate three omnidirectional antennas in a line with 1/2 wavelength between adjancent ones, receiving a signal that arrives at an angle
d = 0.5
Nr = 3
theta_degrees = 20 # direction of arrival
theta = theta_degrees / 180 * np.pi # convert to radians
a = np.exp(-2j * np.pi * d * np.arange(Nr) * np.sin(theta))
print(a)
# we have to do a matrix multiplication of a and tx, so first lets convert both to matrix' instead of numpy arrays which dont let us do 1d matrix math
a = np.asmatrix(a)
tx = np.asmatrix(tx)
# so how do we use this? simple:
r = a.T @ tx # matrix multiply. dont get too caught up by the transpose a, the important thing is we're multiplying the array factor by the tx signal
print(r.shape) # r is now going to be a 2D array, 1d is time and 1d is spatial
# Plot the real part of the first 200 samples of all three elements
fig, (ax1) = plt.subplots(1, 1, figsize=(7, 3))
ax1.plot(np.asarray(r[0, :]).squeeze().real[
0:200]) # the asarray and squeeze are just annoyances we have to do because we came from a matrix
ax1.plot(np.asarray(r[1, :]).squeeze().real[0:200])
ax1.plot(np.asarray(r[2, :]).squeeze().real[0:200])
ax1.set_ylabel("Samples")
ax1.set_xlabel("Time")
ax1.grid()
ax1.legend(['0', '1', '2'], loc=1)
plt.show()
# note the phase shifts, they are also there on the imaginary portions of the samples
# So far this has been simulating the recieving of a signal from a certain angle of arrival
# in your typical DOA problem you are given samples and have to estimate the angle of arrival(s)
# there are also problems where you have multiple receives signals from different directions and one is the SOI while another might be a jammer or interferer you have to null out
# One thing we didnt both doing- lets add noise to this recieved signal.
# AWGN with a phase shift applied is still AWGN so we can add it after or before the array factor is applied, doesnt really matter, we'll do it after
# we need to make sure each element gets an independent noise signal added
n = np.random.randn(Nr, N) + 1j * np.random.randn(Nr, N)
r = r + 0.1 * n
fig, (ax1) = plt.subplots(1, 1, figsize=(7, 3))
ax1.plot(np.asarray(r[0, :]).squeeze().real[
0:200]) # the asarray and squeeze are just annoyances we have to do because we came from a matrix
ax1.plot(np.asarray(r[1, :]).squeeze().real[0:200])
ax1.plot(np.asarray(r[2, :]).squeeze().real[0:200])
ax1.set_ylabel("Samples")
ax1.set_xlabel("Time")
ax1.grid()
ax1.legend(['0', '1', '2'], loc=1)
plt.show()
# conventional beamforming
theta_scan = np.linspace(-1 * np.pi, np.pi, 1000) # 100 different thetas between -180 and +180 degrees
results = []
for theta_i in theta_scan:
# print(theta_i)
w = np.asmatrix(np.exp(-2j * np.pi * d * np.arange(Nr) * np.sin(theta_i))) # look familiar?
r_weighted = np.conj(w) @ r # apply our weights corresponding to the direction theta_i
r_weighted = np.asarray(r_weighted).squeeze() # get it back to a normal 1d numpy array
results.append(np.mean(np.abs(r_weighted) ** 2)) # energy detector
print(theta_scan[np.argmax(results)] * 180 / np.pi) # 19.99999999999998
fig, (ax1) = plt.subplots(1, 1, figsize=(7, 3))
ax1.plot(theta_scan * 180 / np.pi, results) # lets plot angle in degrees
ax1.plot([20], [np.max(results)], 'r.')
ax1.text(-5, np.max(results) + 0.7, '20 degrees')
ax1.set_xlabel("Theta [Degrees]")
ax1.set_ylabel("DOA Metric")
ax1.grid()
plt.show()
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_rgrids([0, 2, 4, 6, 8])
ax.set_rlabel_position(22.5) # Move grid labels away from other labels
plt.show()
