feat: 完成阵列构造, 导向矩阵生成和旋转测试

This commit is contained in:
2025-07-04 15:16:03 +08:00
parent bf8d7fd556
commit 4535f4b329
2 changed files with 540 additions and 48 deletions
+312 -48
View File
@@ -11,6 +11,7 @@ class ArrayType(Enum):
NESTED = 'nested'
COPRIME = 'coprime'
MRA = 'mra'
CYLINDRICAL = "cylindrical"
CUSTOM = 'custom'
class BaseArray(ABC):
@@ -57,12 +58,12 @@ class BaseArray(ABC):
ValueError: If input arrays have mismatched shapes or invalid dimensions.
"""
# position-related
positions = np.array(positions)
if positions.ndim != 2:
raise ValueError(f"`positions` must be a 2D array, got ndim={positions.ndim}.")
# element positions
positions = np.array(positions, dtype=float)
if positions.ndim == 1:
positions = positions.reshape(1, -1) # reshape to row vector
if positions.shape[0] not in (1, 2, 3):
raise ValueError(f"`positions` must be in 1D, 2D, or 3D space. Got dim={positions.shape[0]}.")
raise ValueError(f"positions must be in 1D, 2D, or 3D space. Got dim={positions.shape[0]}.")
self.positions = positions
self.dim = positions.shape[0]
self.count = positions.shape[1]
@@ -75,29 +76,31 @@ class BaseArray(ABC):
self.indices = indices
# errors and perturbations
# TODO: mutual coupling needs to be modeled
# TODO: if position error is provided, the dimension of the array should change w.s.t. the dimension of the error
if position_errors is not None:
position_errors = np.array(position_errors)
position_errors = np.array(position_errors, dtype=float)
if position_errors.shape != positions.shape:
raise ValueError(f"Shape mismatch: position_errors {position_errors.shape} vs positions {positions.shape}")
self.position_errors = position_errors
else:
self.position_errors = np.zeros_like(positions)
self.position_errors = position_errors # None
if gain_errors is not None:
gain_errors = np.array(gain_errors)
gain_errors = np.array(gain_errors, dtype=float)
if gain_errors.shape != positions.shape:
raise ValueError(f"Shape mismatch: gain_errors {gain_errors.shape} vs positions {positions.shape}")
self.gain_errors = gain_errors
else:
self.gain_errors = np.zeros_like(positions)
self.gain_errors = gain_errors # None
if phase_errors is not None:
phase_errors = np.array(phase_errors)
phase_errors = np.array(phase_errors, dtype=float)
if phase_errors.shape != positions.shape:
raise ValueError(f"Shape mismatch: phase_errors {phase_errors.shape} vs positions {positions.shape}")
self.phase_errors = phase_errors
else:
self.phase_errors = np.zeros_like(positions)
self.phase_errors = phase_errors # None
# other attributes
self.spacing = spacing
@@ -223,7 +226,7 @@ class Array1D(BaseArray):
ValueError: If dimension of directions is not 1 or 2.
"""
directions = np.array(directions)
if directions.ndim == 1:
if directions.ndim in (0, 1):
directions = np.reshape(directions, (1, -1))
elif directions.ndim == 2:
if directions.shape[0] != 2:
@@ -251,13 +254,16 @@ class Array1D(BaseArray):
DA = 1j * 2 * np.pi / wavelength * np.outer(positions, np.cos(directions)) * A
# apply gain/phase errors
A *= np.exp(1j * self.phase_errors.flatten())
A *= self.gain_errors.flatten()
if DA is not None:
A_phase = np.exp(1j * self.phase_errors.flatten())[:, None]
if self.gain_errors is not None:
A_gain = self.gain_errors.flatten()[:, None]
DA *= A_phase
DA *= A_gain
A *= A_gain
if DA is not None:
DA *= A_gain
if self.phase_errors is not None:
A_phase = np.exp(1j * self.phase_errors.flatten())[:, None]
A *= A_phase
if DA is not None:
DA *= A_phase
return A, DA
@@ -297,31 +303,36 @@ class Array1D(BaseArray):
ValueError: If `center` is provided but its dimension does not match `target_normal`.
"""
target_normal = np.reshape(np.array(target_normal, dtype=float), (-1, 1))
dim = target_normal.shape[0]
if dim not in (2, 3):
target_normal = (np.array(target_normal, dtype=float)).flatten()
target_dim = target_normal.shape[0]
if target_dim not in (2, 3):
raise ValueError(f"Target normal must be a 2D or 3D vector, got shape {target_normal.shape}.")
positions = self.positions
elem_num = positions.shape[1]
elem_num = self.count
# Determine rotation center
if center is None:
center = np.mean(positions, axis=1)
center = np.reshape(np.array(center, dtype=float), (-1, 1))
if center.shape[0] != dim:
if center.shape[0] != positions.shape[0]:
raise ValueError(
f"Center and target_normal must have the same dimension. "
f"Got center dim = {center.shape[0]}, target_normal dim = {dim}."
f"Array center and element positions must have the same dimension. "
f"Got center dim = {center.shape[0]}, positions dim = {positions.shape[0]}."
)
# Centered positions aligned along y-axis
if dim == 2:
if target_dim == 2:
# Arrange elements along Y-axis and center at origin
centered_pos = np.vstack((np.zeros((1, elem_num)), positions - center))
# Get normal vector
if current_normal is None:
current_normal = np.array([1.0, 0.0]) # assume
else:
current_normal = (np.array(current_normal, dtype=float)).flatten()
current_dim = current_normal.shape[0]
if current_dim != target_dim:
raise ValueError(f"Dimension of target normal ({target_dim}) and current normal ({current_dim}) must match.")
current_normal = current_normal / np.linalg.norm(current_normal)
target_normal = target_normal / np.linalg.norm(target_normal)
# Get the angle between current and target normal vector
@@ -330,7 +341,13 @@ class Array1D(BaseArray):
[np.cos(angle), -np.sin(angle)],
[np.sin(angle), np.cos(angle)]
])
positions_rot = rot_mat @ centered_pos + center[:, np.newaxis]
# expand the dimension of rotation center
center_exp = np.vstack([
np.zeros((1, 1)),
center.reshape(1, 1)
])
# perform rotation and translation
positions_rot = rot_mat @ centered_pos + center_exp
return Array2D(
positions = positions_rot,
type = self.type,
@@ -340,8 +357,8 @@ class Array1D(BaseArray):
gain_errors = self.gain_errors,
phase_errors = self.phase_errors
), rot_mat
elif dim == 3:
# Arrange elements along xOy plane and center at origin
elif target_dim == 3:
# Arrange elements along y-axis and center at origin
centered_pos = np.vstack((
np.zeros((1, elem_num)),
positions - center,
@@ -350,6 +367,11 @@ class Array1D(BaseArray):
# Get normal vector
if current_normal is None:
current_normal = np.array([0, 0, 1])
else:
current_normal = (np.array(current_normal, dtype=float)).flatten()
current_dim = current_normal.shape[0]
if current_dim != target_dim:
raise ValueError(f"Dimension of target normal ({target_dim}) and current normal ({current_dim}) must match.")
current_normal = current_normal / np.linalg.norm(current_normal)
target_normal = target_normal / np.linalg.norm(target_normal)
# Rotate to a 3D direction using axis-angle method
@@ -370,7 +392,15 @@ class Array1D(BaseArray):
[-axis[1], axis[0], 0]
])
rot_mat = np.eye(3) + np.sin(angle) * K + (1 - np.cos(angle)) * K @ K
positions_rot = rot_mat @ centered_pos + center[:, np.newaxis]
# expand the dimension of rotation center
center_exp = np.vstack([
np.zeros((1, 1)),
center.reshape(1, 1),
np.zeros((1, 1))
])
# perform rotation and translation
positions_rot = rot_mat @ centered_pos + center
return Array3D(
positions = positions_rot,
type = self.type,
@@ -403,14 +433,13 @@ class Array2D(BaseArray):
Args:
wavelength (float): Signal wavelength.
directions (np.ndarray): Input directions, must be in one of the following formats:
- shape (D,) or (1, D): D broadside angles (in radians), valid for 1D array
- shape (2, D): D pairs of azimuth and elevation angles (in radians), used for projection to broadside
Returns:
A (np.ndarray): Steering matrix of shape (M, D)
Raises:
ValueError: If dimension of directions is not 1 or 2.
ValueError: If dimension of directions is 2.
"""
directions = np.array(directions)
@@ -422,12 +451,12 @@ class Array2D(BaseArray):
)
az, el = directions[0, :], directions[1, :]
cos_el = np.cos(el)
x_proj = cos_el * np.sin(az)
y_proj = cos_el * np.cos(az)
x_proj = cos_el * np.cos(az)
y_proj = cos_el * np.sin(az)
else:
raise ValueError(
f"Invalid shape for directions ({directions.ndim}). "
"Expected shape (D,) or (1, D) for 1D broadside angles or shape (2, D) for azimuth/elevation."
"Expected shape (2, D) for azimuth/elevation."
)
positions = self.positions
@@ -438,15 +467,115 @@ class Array2D(BaseArray):
A = np.exp(1j * phase)
# apply gain/phase errors
A *= np.exp(1j * self.phase_errors.flatten())
A *= self.gain_errors.flatten()
if self.gain_errors is not None:
A *= self.gain_errors.flatten()
if self.phase_errors is not None:
A *= np.exp(1j * self.phase_errors.flatten())
return A
def rotate_to(self, target_dim: int, rotation_matrix: np.ndarray) -> BaseArray:
assert target_dim == 3 and rotation_matrix.shape == (3, 2)
new_positions = self.element_positions @ rotation_matrix.T
return Array3D(...)
def rotate(
self,
target_normal: np.ndarray,
current_normal: Optional[np.ndarray] = None,
center: Optional[np.ndarray] = None
) -> BaseArray:
"""
Rotate a 2D array to align its normal vector with a target direction.
The array is assumed to be initially aligned along xOy plane (in 3D). The method performs a rigid-body rotation to align the
array's normal with the specified `target_normal`.
The rotation is performed in three steps:
1. Translate the centroid of the array to the origin.
2. Apply the rotation based on the angle between current and target normal vectors.
3. Translate the array back to the specified or original centroid.
Args:
target_normal (np.ndarray): Target direction to align the array's normal to.
Must be a 3D vector of shape (3,).
current_normal (Optional[np.ndarray]): Original normal direction of the array.
Defaults to [0, 0, 1] in 3D.
center (Optional[np.ndarray]): Rotation center, defaults to the array's centroid.
Must be the same dimension as `target_normal`.
Returns:
Tuple[BaseArray, np.ndarray]:
- A new rotated array object (Array3D) with updated element positions.
- The rotation matrix used in the transformation (3x3).
Raises:
ValueError: If `target_normal` is not 3D.
ValueError: If `center` is provided but its dimension does not match `target_normal`.
"""
target_normal = (np.array(target_normal, dtype=float)).flatten()
target_dim = target_normal.shape[0]
if target_dim != 3:
raise ValueError(f"Target normal must be a 3D vector for 2D array, got shape {target_normal.shape}.")
positions = self.positions
elem_num = self.count
# Determine rotation center
if center is None:
center = np.mean(positions, axis=1)
center = np.reshape(np.array(center, dtype=float), (-1, 1))
if center.shape[0] != positions.shape[0]:
raise ValueError(
f"Array center and element positions must have the same dimension. "
f"Got center dim = {center.shape[0]}, positions dim = {positions.shape[0]}."
)
# Arrange elements along xOy plane and center at origin
centered_pos = np.vstack((
positions - center,
np.zeros((1, elem_num))
))
# Get normal vector
if current_normal is None:
current_normal = np.array([0, 0, 1])
else:
current_normal = (np.array(current_normal, dtype=float)).flatten()
current_dim = current_normal.shape[0]
if current_dim != target_dim:
raise ValueError(f"Dimension of target normal ({target_dim}) and current normal ({current_dim}) must match.")
current_normal = current_normal / np.linalg.norm(current_normal)
target_normal = target_normal / np.linalg.norm(target_normal)
# Rotate to a 3D direction using axis-angle method
axis = np.cross(current_normal, target_normal)
norm_axis = np.linalg.norm(axis)
if norm_axis < 1e-8:
# If vectors are parallel or anti-parallel
if np.dot(current_normal, target_normal) > 0:
rot_mat = np.eye(3) # No rotation
else:
rot_mat = np.diag([-1, -1, 1]) # 180-degree rotation
else:
axis = axis / norm_axis
angle = np.arccos(np.dot(current_normal, target_normal))
K = np.array([
[0, -axis[2], axis[1]],
[axis[2], 0, -axis[0]],
[-axis[1], axis[0], 0]
])
rot_mat = np.eye(3) + np.sin(angle) * K + (1 - np.cos(angle)) * K @ K
# expand the dimension of rotation center
center_exp = np.vstack([
center.reshape(2, 1),
np.zeros((1, 1))
])
# perform rotation and translation
positions_rot = rot_mat @ centered_pos + center_exp
return Array3D(
positions = positions_rot,
type = self.type,
spacing = self.spacing,
indices = self.indices,
position_errors = self.position_errors,
gain_errors = self.gain_errors,
phase_errors = self.phase_errors
), rot_mat
class Array3D(BaseArray):
"""
@@ -457,16 +586,151 @@ class Array3D(BaseArray):
super().__init__(*args, **kwargs)
assert self.dim == 3, "Array3D requires 3D element positions."
def steering_matrix(self, wavelength: float, direction: np.ndarray) -> np.ndarray:
"""
Compute 3D steering vector.
def steering_matrix(
self,
wavelength: float,
directions: np.ndarray,
) -> np.ndarray:
"""Compute the steering matrix of 2D array.
Args:
wavelength (float): Wavelength.
direction (np.ndarray): 3D unit vector.
wavelength (float): Signal wavelength.
directions (np.ndarray): Input directions, must be in one of the following formats:
- shape (2, D): D pairs of azimuth and elevation angles (in radians), used for projection to broadside
Returns:
np.ndarray: Complex steering vector.
A (np.ndarray): Steering matrix of shape (M, D)
Raises:
ValueError: If dimension of directions is 2.
"""
directions = np.array(directions)
if directions.ndim == 2:
if directions.shape[0] != 2:
raise ValueError(
f"If directions is 2D, its shape must be (2, D), "
f"but got shape {directions.shape}."
)
az, el = directions[0, :], directions[1, :]
cos_el = np.cos(el)
x_proj = cos_el * np.cos(az)
y_proj = cos_el * np.sin(az)
z_proj = np.sin(el)
else:
raise ValueError(
f"Invalid shape for directions ({directions.ndim}). "
"Expected shape (2, D) for azimuth/elevation."
)
positions = self.positions
# Compute phase shifts
k = 2 * np.pi / wavelength
return np.exp(1j * k * self.element_positions @ direction)
phase = k * (positions[0, :, None] * x_proj + positions[1, :, None] * y_proj
+ positions[2, :, None] * z_proj)
A = np.exp(1j * phase)
# apply gain/phase errors
if self.gain_errors is not None:
A *= self.gain_errors.flatten()
if self.phase_errors is not None:
A *= np.exp(1j * self.phase_errors.flatten())
return A
def rotate(
self,
target_normal: np.ndarray,
current_normal: Optional[np.ndarray] = None,
center: Optional[np.ndarray] = None
) -> BaseArray:
"""
Rotate a 3D array to align its normal vector with a target direction.
The method performs a rigid-body rotation to align the array's normal with the specified `target_normal`.
The rotation is performed in three steps:
1. Translate the centroid of the array to the origin.
2. Apply the rotation based on the angle between current and target normal vectors.
3. Translate the array back to the specified or original centroid.
Args:
target_normal (np.ndarray): Target direction to align the array's normal to.
Must be a 3D vector of shape (3,).
current_normal (Optional[np.ndarray]): Original normal direction of the array.
Defaults to [0, 0, 1] in 3D.
center (Optional[np.ndarray]): Rotation center, defaults to the array's centroid.
Must be the same dimension as `target_normal`.
Returns:
Tuple[BaseArray, np.ndarray]:
- A new rotated array object (Array3D) with updated element positions.
- The rotation matrix used in the transformation (3x3).
Raises:
ValueError: If `target_normal` is not 3D.
ValueError: If `center` is provided but its dimension does not match element `positions`.
"""
target_normal = (np.array(target_normal, dtype=float)).flatten()
target_dim = target_normal.shape[0]
if target_dim != 3:
raise ValueError(f"Target normal must be a 3D vector for 2D array, got shape {target_normal.shape}.")
positions = self.positions
elem_num = self.count
# Determine rotation center
if center is None:
center = np.mean(positions, axis=1)
center = np.reshape(np.array(center, dtype=float), (-1, 1))
if center.shape[0] != positions.shape[0]:
raise ValueError(
f"Array center and element positions must have the same dimension. "
f"Got center dim = {center.shape[0]}, positions dim = {positions.shape[0]}."
)
# Arrange elements along xOy plane and center at origin
centered_pos = positions - center
# Get normal vector
if current_normal is None:
current_normal = np.array([0, 0, 1])
else:
current_normal = (np.array(current_normal, dtype=float)).flatten()
current_dim = current_normal.shape[0]
if current_dim != target_dim:
raise ValueError(f"Dimension of target normal ({target_dim}) and current normal ({current_dim}) must match.")
current_normal = current_normal / np.linalg.norm(current_normal)
target_normal = target_normal / np.linalg.norm(target_normal)
# Rotate to a 3D direction using axis-angle method
axis = np.cross(current_normal, target_normal)
norm_axis = np.linalg.norm(axis)
if norm_axis < 1e-8:
# If vectors are parallel or anti-parallel
if np.dot(current_normal, target_normal) > 0:
rot_mat = np.eye(3) # No rotation
else:
rot_mat = np.diag([-1, -1, 1]) # 180-degree rotation
else:
axis = axis / norm_axis
angle = np.arccos(np.dot(current_normal, target_normal))
K = np.array([
[0, -axis[2], axis[1]],
[axis[2], 0, -axis[0]],
[-axis[1], axis[0], 0]
])
rot_mat = np.eye(3) + np.sin(angle) * K + (1 - np.cos(angle)) * K @ K
# expand the dimension of rotation center
center_exp = center.reshape(3, 1)
# perform rotation and translation
positions_rot = rot_mat @ centered_pos + center_exp
return Array3D(
positions = positions_rot,
type = self.type,
spacing = self.spacing,
indices = self.indices,
position_errors = self.position_errors,
gain_errors = self.gain_errors,
phase_errors = self.phase_errors
), rot_mat