Skip to content

Commit

Permalink
all comments in English now
Browse files Browse the repository at this point in the history
  • Loading branch information
zhiim committed Jan 20, 2024
1 parent d88f49c commit cc2ca28
Show file tree
Hide file tree
Showing 11 changed files with 174 additions and 158 deletions.
32 changes: 18 additions & 14 deletions classical_doa/algorithm/broadband/cssm.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,14 +10,15 @@ def cssm(received_data, num_signal, array_position, fs, angle_grids, fre_ref,
"""Coherent Signal Subspace Method (CSSM) for wideband DOA estimation.
Args:
received_data : 阵列接受信号
num_signal : 信号个数
array_position : 阵元位置, 应该是numpy array的形式, 行向量列向量均可
fs: 采样频率
angle_grids : 空间谱的网格点, 应该是numpy array的形式
fre_ref: 参考频点
pre_estimate: 角度预估计值
unit : 角度的单位制, `rad`代表弧度制, `deg`代表角度制. Defaults to
received_data : Array received signals
num_signal : Number of signals
array_position : Position of array elements. It should be a numpy array
fs: sampling frequency
angle_grids : Angle grids corresponding to spatial spectrum. It should
be a numpy array.
fre_ref: reference frequency
pre_estimate: pre-estimated angles
unit : Unit of angle, 'rad' for radians, 'deg' for degrees. Defaults to
'deg'.
References:
Expand All @@ -34,24 +35,27 @@ def cssm(received_data, num_signal, array_position, fs, angle_grids, fre_ref,
pre_estimate = pre_estimate.reshape(1, -1)
array_position = array_position.reshape(-1, 1)

# 频域下的阵列接收信号
# Divide the received signal into multiple frequency points
signal_fre_bins = np.fft.fft(received_data, axis=1)
fre_bins = np.fft.fftfreq(num_snapshots, 1 / fs)

# 计算参考频点下,预估计角度对应的流型矩阵
# Calculate the manifold matrix corresponding to the pre-estimated angles at
# the reference frequency point
matrix_a_ref = np.exp(-1j * 2 * np.pi * fre_ref / C *\
array_position @ pre_estimate)

for i, fre in enumerate(fre_bins):
# 每个频点下,角度预估值对应的流型矩阵
# Manifold matrix corresponding to the pre-estimated angles at
# each frequency point
matrix_a_f = np.exp(-1j * 2 * np.pi * fre / C *\
array_position @ np.sin(pre_estimate))
matrix_q = matrix_a_f @ matrix_a_ref.transpose().conj()
# 对matrix_q进行奇异值分解
# Perform singular value decomposition on matrix_q
matrix_u, _, matrix_vh = np.linalg.svd(matrix_q)
# RSS法构造最佳聚焦矩阵
# Construct the optimal focusing matrix using the RSS method
matrix_t_f = matrix_vh.transpose().conj() @ matrix_u.transpose().conj()
# 将每个频点对应的接收信号聚焦到参考频点
# Focus the received signals at each frequency point to the reference
# frequency point
signal_fre_bins[:, i] = matrix_t_f @ signal_fre_bins[:, i]

spectrum = music(received_data=signal_fre_bins, num_signal=num_signal,
Expand Down
19 changes: 10 additions & 9 deletions classical_doa/algorithm/broadband/imusic.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,13 +9,15 @@ def imusic(received_data, num_signal, array_position, fs, angle_grids,
"""Incoherent MUSIC estimator for wideband DOA estimation.
Args:
received_data : 阵列接受信号
num_signal : 信号个数
array_position : 阵元位置, 应该是numpy array的形式, 行向量列向量均可
fs: 采样频率
angle_grids : 空间谱的网格点, 应该是numpy array的形式
num_groups: FFT的组数, 每一组都独立做FFT
unit : 角度的单位制, `rad`代表弧度制, `deg`代表角度制. Defaults to
received_data : Array received signals
num_signal : Number of signals
array_position : Position of array elements. It should be a numpy array
fs: sampling frequency
angle_grids : Angle grids corresponding to spatial spectrum. It should
be a numpy array.
num_groups: Divide sampling points into serveral groups, and do FFT
separately in each group
unit : Unit of angle, 'rad' for radians, 'deg' for degrees. Defaults to
'deg'.
References:
Expand All @@ -27,7 +29,7 @@ def imusic(received_data, num_signal, array_position, fs, angle_grids,
signal_fre_bins, fre_bins = divide_into_fre_bins(received_data, num_groups,
fs)

# 对每一个频点运行MUSIC算法
# MUSIC algorithm in every frequency point
spectrum_fre_bins = np.zeros((signal_fre_bins.shape[1], angle_grids.size))
for i, fre in enumerate(fre_bins):
spectrum_fre_bins[i, :] = music(received_data=signal_fre_bins[:, i, :],
Expand All @@ -37,7 +39,6 @@ def imusic(received_data, num_signal, array_position, fs, angle_grids,
angle_grids=angle_grids,
unit=unit)

# 取平均
spectrum = np.mean(spectrum_fre_bins, axis=0)

return np.squeeze(spectrum)
34 changes: 18 additions & 16 deletions classical_doa/algorithm/broadband/tops.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,15 +15,16 @@ def tops(received_data, num_signal, array_position, fs, num_groups, angle_grids,
DOA estimation.
Args:
received_data : 阵列接受信号
num_signal : 信号个数
array_position : 阵元位置, 应该是numpy array的形式, 行向量列向量均可
fs: 采样频率
num_groups: FFT的组数, 每一组都独立做FFT
angle_grids : 空间谱的网格点, 应该是numpy array的形式
fre_ref: 参考频点
unit : 角度的单位制, `rad`代表弧度制, `deg`代表角度制. Defaults to
'deg'.
received_data: received signals from the array.
num_signal: Number of signals.
array_position: Array element positions, should be a numpy array
fs: Sampling frequency.
num_groups: Number of groups for FFT, each group performs an
independent FFT.
angle_grids: Grid points of spatial spectrum, should be a numpy array.
fre_ref: Reference frequency point.
unit: Unit of angle measurement, 'rad' for radians, 'deg' for degrees.
Defaults to 'deg'.
References:
Yoon, Yeo-Sun, L.M. Kaplan, and J.H. McClellan. “TOPS: New DOA Estimator
Expand All @@ -39,7 +40,7 @@ def tops(received_data, num_signal, array_position, fs, num_groups, angle_grids,
signal_fre_bins, fre_bins = divide_into_fre_bins(received_data, num_groups,
fs)

# index of reference frequency if FFT output
# index of reference frequency in FFT output
ref_index = int(fre_ref / (fs / fre_bins.size))
# get signal space of reference frequency
signal_space_ref = get_signal_space(signal_fre_bins[:, ref_index, :],
Expand All @@ -50,33 +51,34 @@ def tops(received_data, num_signal, array_position, fs, num_groups, angle_grids,
matrix_d = np.empty((num_signal, 0), dtype=np.complex_)

for j, fre in enumerate(fre_bins):
# 计算当前频点对应的噪声子空间
# calculate noise subspace for the current frequency point
noise_space_f = get_noise_space(signal_fre_bins[:, j, :],
num_signal)

# 构造变换矩阵
# construct transformation matrix
matrix_phi = np.exp(-1j * 2 * np.pi * (fre - fre_ref) / C *
array_position * np.sin(grid))
matrix_phi = np.diag(np.squeeze(matrix_phi))

# 使用变换矩阵将参考频点的信号子空间变换到当前频点
# transform the signal subspace of the reference frequency to the
# current frequency using the transformation matrix
matrix_u = matrix_phi @ signal_space_ref

# 构造投影矩阵,减小矩阵U中的误差
# construct projection matrix to reduce errors in matrix U
matrix_a_f = np.exp(-1j * 2 * np.pi * fre / C *
array_position * np.sin(grid))
matrix_p = np.eye(num_antennas) -\
1 / (matrix_a_f.transpose().conj() @ matrix_a_f) *\
matrix_a_f @ matrix_a_f.transpose().conj()

# 使用投影矩阵对矩阵U进行投影
# project matrix U using the projection matrix
matrix_u = matrix_p @ matrix_u

matrix_d = np.concatenate((matrix_d,
matrix_u.T.conj() @ noise_space_f),
axis=1)

# 使用矩阵D中的最小特征值构造空间谱
# construct spatial spectrum using the minimum eigenvalue of matrix D
_, s, _ = np.linalg.svd(matrix_d)
spectrum[i] = 1 / min(s)

Expand Down
17 changes: 9 additions & 8 deletions classical_doa/algorithm/esprit.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,12 +10,12 @@ def esprit(received_data, num_signal, array_position, signal_fre, unit="deg"):
the reference paper.
Args:
received_data : 阵列接受信号
num_signal : 信号个数
array_position : 阵元位置, 应该是numpy array的形式, 行向量列向量均可
signal_fre: 信号频率
unit : 返回的估计角度的单位制, `rad`代表弧度制, `deg`代表角度制.
Defaults to 'deg'.
received_data : Array received signals
num_signal : Number of signals
array_position : Position of array elements. It should be a numpy array
signal_fre: Signal frequency
unit : Unit of angle, 'rad' for radians, 'deg' for degrees. Defaults to
'deg'.
Reference:
Roy, R., and T. Kailath. “ESPRIT-Estimation of Signal Parameters via
Expand All @@ -28,15 +28,16 @@ def esprit(received_data, num_signal, array_position, signal_fre, unit="deg"):
# get signal space of two sub array. Each sub array consists of M-1 antennas
matrix_e_x = signal_space[:-1, :]
matrix_e_y = signal_space[1:, :]
# 两个子阵列中对应阵元之间的固定间距确保了旋转不变性
# the fixed distance of corresponding elements in two sub-array ensures
# the rotational invariance
sub_array_spacing = array_position[1] - array_position[0]

matrix_c = np.hstack((matrix_e_x, matrix_e_y)).transpose().conj() @\
np.hstack((matrix_e_x, matrix_e_y))

# get eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(matrix_c)
sorted_index = np.argsort(np.abs(eigenvalues))[::-1] # 由大到小排序的索引
sorted_index = np.argsort(np.abs(eigenvalues))[::-1] # descending order
matrix_e = eigenvectors[:, sorted_index[:2 * num_signal]]

# take the upper right and lower right sub matrix
Expand Down
49 changes: 29 additions & 20 deletions classical_doa/algorithm/music.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,45 +10,49 @@ def music(received_data, num_signal, array_position, signal_fre, angle_grids,
"""1D MUSIC
Args:
received_data : 阵列接受信号
num_signal : 信号个数
array_position : 阵元位置, 应该是numpy array的形式, 行向量列向量均可
signal_fre: 信号频率
angle_grids : 空间谱的网格点, 应该是numpy array的形式
unit : 角度的单位制, `rad`代表弧度制, `deg`代表角度制. Defaults to
received_data : Array received signals
num_signal : Number of signals
array_position : Position of array elements. It should be a numpy array
signal_fre: Signal frequency
angle_grids : Angle grids corresponding to spatial spectrum. It should
be a numpy array.
unit : Unit of angle, 'rad' for radians, 'deg' for degrees. Defaults to
'deg'.
"""
noise_space = get_noise_space(received_data, num_signal)

# 变为列向量用于后续矩阵计算
# Reshape to column vector for matrix calculations
array_position = np.reshape(array_position, (-1, 1))

if unit == "deg":
angle_grids = angle_grids / 180 * np.pi
angle_grids = np.reshape(angle_grids, (1, -1))

# 计算所有网格点信号入射时的流型矩阵
# Calculate the manifold matrix when there are incident signal in all
# grid points
tau_all_grids = 1 / C * array_position @ np.sin(angle_grids)
manifold_all_grids = np.exp(-1j * 2 * np.pi * signal_fre * tau_all_grids)

v = noise_space.transpose().conj() @ manifold_all_grids

# 矩阵v的每一列对应一个入射信号, 对每一列求二范数的平方
# Each column of matrix v corresponds to an incident signal, calculate the
# square of the 2-norm for each column
spectrum = 1 / np.linalg.norm(v, axis=0) ** 2

return np.squeeze(spectrum)


def root_music(received_data, num_signal, array_position, signal_fre,
unit="deg"):
"""Root-MUSIC
Args:
received_data : 阵列接受信号
num_signal : 信号个数
array_position : 阵元位置, 应该是numpy array的形式, 行向量列向量均可
signal_fre: 信号频率
unit : 角度的单位制, `rad`代表弧度制, `deg`代表角度制. Defaults to
'deg'.
received_data : Array of received signals
num_signal : Number of signals
array_position : Position of array elements. It should be a numpy array
signal_fre: Signal frequency
unit: The unit of the angle, `rad` represents radian, `deg` represents
degree. Defaults to 'deg'.
References:
Rao, B.D., and K.V.S. Hari. “Performance Analysis of Root-Music.”
Expand All @@ -60,9 +64,14 @@ def root_music(received_data, num_signal, array_position, signal_fre,
num_antennas = array_position.size
antenna_spacing = array_position[1] - array_position[0]

# 由于numpy提供的多项式求解函数需要以多项式的系数作为输入,而系数的提取非常
# 复杂, 此处直接搬用doatools中rootMMUSIC的实现代码
# 也可以使用sympy库求解多项式方程,但是计算量会更大
# Since the polynomial solving function provided by numpy requires the
# coefficients of the polynomial as input, and extracting the coefficients
# is very complex, so the implementation code of rootMMUSIC in doatools is
# directly used here.

# Alternatively, the sympy library can be used to solve polynomial
# equations, but it will be more computationally expensive.

# Compute the coefficients for the polynomial.
matrix_c = noise_space @ noise_space.transpose().conj()
coeff = np.zeros((num_antennas - 1,), dtype=np.complex_)
Expand All @@ -72,8 +81,8 @@ def root_music(received_data, num_signal, array_position, signal_fre,
# Find the roots of the polynomial.
z = np.roots(coeff)

# 为了防止同时取到一对共轭根, 只取单位圆内的根
# find k roots inside and closest to the unit circle
# To avoid simultaneously obtaining a pair of complex conjugate roots, only
# take roots inside the unit circle
roots_inside_unit_circle = np.extract(np.abs(z) <= 1, z)
sorted_index = np.argsort(np.abs(np.abs(roots_inside_unit_circle) - 1))
chosen_roots = roots_inside_unit_circle[sorted_index[:num_signal]]
Expand Down
13 changes: 7 additions & 6 deletions classical_doa/algorithm/sparse.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,12 +8,13 @@ def omp(received_data, num_signal, array_position, signal_fre, angle_grids,
"""OMP based sparse representation algorithms for DOA estimation
Args:
received_data : 阵列接受信号
num_signal : 信号个数
array_position : 阵元位置, 应该是numpy array的形式, 行向量列向量均可
signal_fre: 信号频率
angle_grids : 空间谱的网格点, 应该是numpy array的形式
unit : 角度的单位制, `rad`代表弧度制, `deg`代表角度制. Defaults to
received_data : Array received signals
num_signal : Number of signals
array_position : Position of array elements. It should be a numpy array
signal_fre: Signal frequency
angle_grids : Angle grids corresponding to spatial spectrum. It should
be a numpy array.
unit : Unit of angle, 'rad' for radians, 'deg' for degrees. Defaults to
'deg'.
Reference:
Expand Down
17 changes: 9 additions & 8 deletions classical_doa/algorithm/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ def get_noise_space(received_data, num_signal):
(received_data @ received_data.transpose().conj())

eigenvalues, eigenvectors = np.linalg.eig(corvariance_matrix)
sorted_index = np.argsort(np.abs(eigenvalues)) # 由小到大排序的索引
sorted_index = np.argsort(np.abs(eigenvalues)) # ascending order
noise_space = eigenvectors[:, sorted_index[:-num_signal]]

return noise_space
Expand All @@ -23,7 +23,7 @@ def get_signal_space(received_data, num_signal):
(received_data @ received_data.transpose().conj())

eigenvalues, eigenvectors = np.linalg.eig(corvariance_matrix)
sorted_index = np.argsort(np.abs(eigenvalues)) # 由小到大排序的索引
sorted_index = np.argsort(np.abs(eigenvalues)) # ascending order
noise_space = eigenvectors[:, sorted_index[-num_signal:]]

return noise_space
Expand All @@ -39,21 +39,22 @@ def divide_into_fre_bins(received_data, num_groups, fs):
fs : sampling frequency
Returns:
`signal_fre_bins`: 一个m*n*l维的矩阵, 其中m为阵元数, n为FFT的点数,
l为组数
`fre_bins`: 每一个FFT点对应的频点
`signal_fre_bins`: a m*n*l tensor, in which m equals to number of
antennas, n is equals to point of FFT, l is the number of groups
`fre_bins`: corresponding freqeuncy of each point in FFT output
"""
num_snapshots = received_data.shape[1]

n_each_group = num_snapshots // num_groups # 每一组包含的采样点数
# number of sampling points in each group
n_each_group = num_snapshots // num_groups
if n_each_group < 128:
n_fft = 128 # 如果点数太少, 做FFT时补零
n_fft = 128 # zero padding when sampling points is not enough
else:
n_fft = n_each_group

signal_fre_bins = np.zeros((received_data.shape[0], n_fft, num_groups),
dtype=np.complex_)
# 每一组独立做FFT
# do FTT separately in each group
for group_i in range(num_groups):
signal_fre_bins[:, :, group_i] = np.fft.fft(
received_data[:, group_i * n_each_group:(group_i+1) * n_each_group],
Expand Down
Loading

0 comments on commit cc2ca28

Please sign in to comment.