import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 设置中文字体和负号显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 物理常数
c = 3e8 # 光速 (m/s)
def design_fda(N, M, f0, delta_f, r0, theta0, d=None):
"""
设计频控阵并计算波束汇聚所需的相位补偿
"""
# 1. 计算默认阵元间距 (半波长)
lambda0 = c / f0
d = d if d is not None else lambda0 / 2
# 2. 生成阵元位置 [-Nd, ..., 0, ..., Nd]
n_indices = np.arange(-N, N + 1)
positions = n_indices * d
# 3. 生成频率矩阵 (2N+1) x (M+1)
m_range = np.arange(M + 1)
n_abs = np.abs(n_indices)
# Δf_n = Δf * log(|n| + 1)
delta_f_n = delta_f * np.log(n_abs + 1 + 1e-12) # +1e-12 避免 log(0)
# Δf_m = Δf * log(m+1)
delta_f_m = delta_f * np.log(m_range + 1)
# f_{n,m} = f0 + Δf_m + Δf_n
freq_matrix = f0 + delta_f_n[:, np.newaxis] + delta_f_m[np.newaxis, :]
# 4. 计算目标点直角坐标
target_x = r0 * np.cos(theta0)
target_y = r0 * np.sin(theta0)
# 5. 计算距离和传播时延
distances = np.sqrt((positions - target_x) ** 2 + target_y ** 2)
time_delays = distances / c
# 6. 计算相位补偿矩阵 (补偿传播相位)
phase_comp = 2 * np.pi * freq_matrix * time_delays[:, np.newaxis]
# 7. 添加额外补偿确保目标点同相
ref_phase = phase_comp[N, 0] # 中心阵元第一个载波的相位
phase_comp -= ref_phase # 使参考相位为零
return positions, freq_matrix, phase_comp
def plot_beam_patterns(N, M, f0, delta_f, r0, theta0):
"""
绘制频控阵的波束方向图（含二维俯视图）
"""
# 设计频控阵
positions, freq_matrix, phase_comp = design_fda(N, M, f0, delta_f, r0, theta0)
# 创建计算网格
r_values = np.linspace(0.5 * r0, 2 * r0, 200) # 距离范围
theta_values = np.linspace(-np.pi / 2, np.pi / 2, 200) # 角度范围：-90°到90°
R, Theta = np.meshgrid(r_values, theta_values)
# 转换为直角坐标用于场强计算
X = R * np.cos(Theta)
Y = R * np.sin(Theta)
# 计算网格上每个点的场强
print("开始计算场强分布...")
field_intensity = np.zeros_like(R, dtype=complex)
num_elements = len(positions)
num_freqs = freq_matrix.shape[1]
wavelengths = c / freq_matrix # 预计算所有波长
# 优化计算：向量化操作
for i in range(num_elements):
elem_x = positions[i]
elem_y = 0 # 假设阵列沿x轴放置
# 预计算距离网格
dist_grid = np.sqrt((elem_x - X) ** 2 + (elem_y - Y) ** 2)
for m in range(num_freqs):
# 计算传播相位
prop_phase = 2 * np.pi * dist_grid / wavelengths[i, m]
# 计算场强贡献
field_intensity += np.exp(1j * (phase_comp[i, m] - prop_phase))
# 取场强的模值
beam_pattern = np.abs(field_intensity)
print("场强计算完成！")
# 将角度转换为度
theta_deg = np.degrees(theta_values)
r_km = r_values / 1000
# ================== 2D 距离和角度切面图 ==================
plt.figure(figsize=(14, 6))
# 角度维切面 (固定距离r0)
plt.subplot(1, 2, 1)
r_idx = np.argmin(np.abs(r_values - r0))
plt.plot(theta_deg, beam_pattern[:, r_idx])
plt.axvline(np.degrees(theta0), color='r', linestyle='--', label='目标角度')
plt.xlabel('角度 (度)')
plt.ylabel('场强幅度')
plt.title(f'固定距离 r={r0 / 1000:.1f}km 的角度维方向图')
plt.legend()
plt.grid(True)
plt.xlim(-90, 90)
# 距离维切面 (固定角度theta0)
plt.subplot(1, 2, 2)
angle_idx = np.argmin(np.abs(theta_values - theta0))
plt.plot(r_km, beam_pattern[angle_idx, :])
plt.axvline(r0 / 1000, color='r', linestyle='--', label='目标距离')
plt.xlabel('距离 (km)')
plt.ylabel('场强幅度')
plt.title(f'固定角度 θ={np.degrees(theta0):.1f}° 的距离维方向图')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
# ================== 3D 场强距离角度图 ==================
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# 创建角度和距离网格
THETA, R_GRID = np.meshgrid(theta_deg, r_km)
# 绘制3D曲面图
surf = ax.plot_surface(THETA, R_GRID, beam_pattern.T,
cmap='viridis',
edgecolor='none',
alpha=0.8,
antialiased=True)
# 标记目标点
ax.scatter([np.degrees(theta0)], [r0 / 1000], [np.max(beam_pattern)],
s=50, c='red', marker='o', label='目标点')
ax.set_xlabel('角度 (度)')
ax.set_ylabel('距离 (km)')
ax.set_zlabel('场强幅度')
ax.set_title('频控阵三维距离-角度场强分布')
ax.legend()
# 添加颜色条
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, label='场强幅度')
# 设置视角
ax.view_init(elev=30, azim=45)
plt.tight_layout()
plt.show()
# ================== 2D 俯视图 ==================
plt.figure(figsize=(10, 8))
# 创建俯视图：角度-距离平面上的场强分布
plt.pcolormesh(theta_deg, r_km, beam_pattern.T, shading='gouraud', cmap='viridis')
# 添加轮廓线
levels = np.linspace(0.2 * np.max(beam_pattern), np.max(beam_pattern), 10)
plt.contour(theta_deg, r_km, beam_pattern.T, levels=levels, colors='white', linewidths=0.5, alpha=0.7)
# 设置坐标轴和标签
plt.xlabel('角度 (度)')
plt.ylabel('距离 (km)')
plt.title('频控阵二维俯视场强分布')
plt.colorbar(label='场强幅度')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
# 设置坐标范围
plt.xlim(-90, 90)
plt.ylim(0.5 * r0 / 1000, 2 * r0 / 1000)
plt.tight_layout()
plt.show()
# 运行可视化
if __name__ == "__main__":
# 参数设置优化
N = 5
M = 3
f0 = 2.4e9 # 10 GHz 中心频率
delta_f = 300e3 # 300 kHz 频率步进
# 目标点参数
r0 = 4000 # 1 km 距离
theta0 = np.pi / 4 # 30度方位角
print("开始计算波束方向图...")
# 绘制方向图
plot_beam_patterns(N, M, f0, delta_f, r0, theta0)
print("计算完成！")在这个代码的基础上，修改二维场强图的场强颜色对应柱，现在是一片空白
最新发布
扫一扫
527
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
