# import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import hamming
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei'] # 使用黑体
plt.rcParams['axes.unicode_minus'] = False # 处理负号的显示
# 常量定义
C = 3e8 # 光速
class PULSE_RADAR:
def __init__(self, bw, pa, theta0, f0, pri, T, fs, DR, SNR, RCS,
barker_code=np.array([1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1]),
FM_beta=5, M=10, costas_code=np.array([2, 4, 8, 5, 10, 9, 7, 3, 6, 1]),):
""" 信号参数 """
# 光速
self.c = 3e8
# 脉冲宽度
self.PW = DR*pri
# 信号带宽
self.Bw = bw
# 脉冲幅度
self.PA = pa
# 信号初始相位
self.theta0 = theta0
# 载频
self.f0 = f0
# 波长
self.lambdas = self.c / self.f0
# 脉冲重复周期
self.pri = pri
# 调频斜率
self.K = bw/self.PW
# 回波窗持续时间
self.T = T
# 采样频率
self.fs = fs
# SNR
self.SNR = SNR
# RCS
self.RCS = RCS
# 巴克码序列
self.barker_code = barker_code
# 调频指数
self.FM_beta = FM_beta
# Frank码阶数
self.M = M
# 回波窗内信号采样点数
self.N_effective = int(np.around(fs * self.PW))
# NLFM插值因子
self.int_num = int(self.fs / (2 * self.Bw))
self.interp_coe = 4
# costas跳频序列
self.costas_code = costas_code
def interpulse_modulation(self, modulation_type):
"""
脉间调制
:param modulation_type:
:return:
"""
# if modulation_type == "jagged":
def costas_signal_generate(self):
"""costas基带信号生成"""
N = int(np.around(self.fs * self.PW))
t = np.linspace(0, self.PW, N, endpoint=False)
num_costas = len(self.costas_code)
# 生成costas跳频频点序列
costas_f_list = - self.Bw / 2 + (self.costas_code - 1) * self.Bw / (num_costas - 1)
# 每个频率编码码片的持续时间
T_costaschip = self.PW / num_costas
# 每个频率编码码片持续时间对应的采样点
num_points_per_segment = int(T_costaschip * self.fs)
# 用于生成单个频率编码码片基带信号的时间轴
t_costaschip = np.linspace(0, T_costaschip, num_points_per_segment, endpoint=False)
xt = np.zeros_like(t, dtype=complex)
for i in range(num_costas):
start_index = i * num_points_per_segment
end_index = (i + 1) * num_points_per_segment
xt[start_index:end_index] = np.exp(1j * 2 * np.pi * costas_f_list[i] * t_costaschip)
return xt
def frank_generate(self):
"""
Frank码IQ调制复信号生成生成
"""
N = self.M ** 2
frank_phase = list(None for _ in range(N))
for i in range(self.M):
for j in range(self.M):
frank_phase[j + i * self.M] = 2 * np.pi * j * i / self.M
frank_complex = np.exp(1j * np.array(frank_phase))
return frank_complex
def frank_code_generate(self, td, n):
"""生成Frank码基带信号"""
# 码长
N = self.M ** 2
# 码元宽度
tb = self.PW / N
N_chip = int(np.around(tb * self.fs))
frank_code = np.zeros(N * N_chip, np.complex64)
frank_z = self.frank_generate()
for i in range(N):
frank_code[(i * N_chip):(i * N_chip + N_chip)] = frank_z[i] + frank_code[
(i * N_chip):(i * N_chip + N_chip)]
bool_t = (np.abs(td[n, :]) < self.PW / 2)
coords = np.where(bool_t)
start_index = coords[0][0]
end_index = start_index + frank_code.shape[0]
baseband_signal = np.zeros_like(bool_t, dtype=np.complex64)
baseband_signal[start_index:end_index] = baseband_signal[start_index:end_index] + frank_code
return baseband_signal
def code_to_signal(self, td, n):
"""将巴克码转换为对应的基带信号"""
bits = 2 * self.barker_code - 1
bool_t = (np.abs(td[n, :]) < self.PW / 2)
coords = np.where(bool_t)
start_index = coords[0][0]
end_index = start_index + self.N_effective
code_expand = code_sampling(bits, self.N_effective, self.Bw, self.fs)
baseband_signal = bool_t.astype(int)
baseband_signal[start_index:end_index] = baseband_signal[start_index:end_index] * code_expand
return baseband_signal
def hamming_nlfm_generate(self):
""" 利用hamming窗生成NLFM信号 """
# 计算时间采样间隔和点数
ts = 1 / (2 * self.Bw)
N = int(np.around(self.PW / ts))
# 生成时间向量 t（从0到T-ts，共N点）
t = np.linspace(0, self.PW, N, endpoint=False) # endpoint=False 排除T点，类似MATLAB的0:ts:T-ts
t1 = t - self.PW / 2 # 中心在0的时间向量
# 计算采样频率和频率分辨率
fs = 1 / ts # 采样频率：1/ts = B = 8e6 Hz
fbin = fs / N # 频率分辨率：fs/N
f = np.arange(0, fs, fbin) # 频率向量从0到fs-fbin，共N点
f1 = f - fs / 2 # 中心在0的频率向量
# 生成汉明窗
wa = hamming(N) # 使用Scipy的hamming函数，返回N点汉明窗
# 插值窗函数、频率向量和时间向量到更高密度
# 使用numpy.interp进行线性插值：新x坐标为原索引的插值版本
x_old = np.arange(len(wa)) # 原索引：0到N-1
x_new = np.linspace(0, len(wa) - 1, len(wa) * self.int_num) # 新索引：0到N-1，插值到int_num倍密度
wa_interp = np.interp(x_new, x_old, wa) # 插值窗函数
f1_interp = np.interp(x_new, x_old, f1) # 插值频率向量
t1_interp = np.interp(x_new, x_old, t1) # 插值时间向量
# 计算群延迟 Tg(f)
K = self.PW / np.sum(wa_interp) # 归一化常数
cum_sum_wa = np.cumsum(wa_interp) # 累积和窗函数
Tg = K * cum_sum_wa - self.PW / 2 # 群延迟：Tg(f) = K * 累积和(窗) - T/2
# 进一步插值时间向量 t1_interp 到 interp_coe 倍密度
x_old2 = np.arange(len(t1_interp)) # t1_interp 的索引
x_new2 = np.linspace(0, len(t1_interp) - 1, len(t1_interp) * self.interp_coe) # 新索引
t2_interp = np.interp(x_new2, x_old2, t1_interp) # 插值时间向量
# 计算频率函数 f(t)：通过插值群延迟的逆函数得到
# 使用 np.interp：在 Tg 上插值 f1_interp，求值点为 t2_interp（假设 Tg 单调）
fq1 = np.interp(t2_interp, Tg, f1_interp) # fq1 是插值后的频率序列
# 降采样 fq1：从索引1开始（Python索引0-based），每隔 interp_coe 取一个点
fq = fq1[1::self.interp_coe] # 类似MATLAB的 fq1(2:interp_coe:end)
# 计算相位：相位是频率的积分，近似为累积和
dt = ts / self.int_num # 时间步长：原始 ts 除以插值倍数
phase_w = 2 * np.pi * np.cumsum(fq) * dt # 相位 = 2π * ∫f(t)dt ≈ 2π * sum(fq) * dt
xt = np.exp(1j * phase_w) # 信号：e^(j*phase_w)
return xt
def signal_generate(self, pulse_num, tr, td, echo_signal, R, Rmin, Rmax, Nwind, v, signal_type):
""" 回波信号生成 """
for n in range(pulse_num):
tr[n, :] = np.linspace(2 * Rmin / C, 2 * Rmax / C, Nwind) # 可检测到的时间范围（由距离决定），时间轴
td[n, :] = tr[n, :] - 2 * (R - v * n * self.pri) / C # pulse_num 个 PRT 的回波时间变量，分别存储在矩阵 td 的每一行
if signal_type == 1:
echo_signal[n, :] = self.RCS * self.PA * (np.abs(td[n, :]) < self.PW / 2) * np.exp(
1j * 2 * np.pi * self.f0 * td[n, :] + 1j * np.pi * self.K * td[n, :] ** 2 + 1j * self.theta0)
elif signal_type == 2:
baseband_signal = self.code_to_signal(td, n)
# echo_signal[n, :] = (self.RCS * self.PA * (np.abs(td[n, :]) < self.PW / 2) * baseband_signal *
# np.exp(1j * 2 * np.pi * self.f0 * td[n, :] + 1j * self.theta0))
echo_signal[n, :] = (self.RCS * self.PA * baseband_signal * np.exp(1j * 2 * np.pi * self.f0 * td[n, :] + 1j * self.theta0))
elif signal_type == 3:
"""Frank码信号生成"""
baseband_signal = self.frank_code_generate(td, n)
echo_signal[n, :] = (self.RCS * self.PA * baseband_signal * np.exp(
1j * 2 * np.pi * self.f0 * td[n, :] + 1j * self.theta0))
elif signal_type == 4:
"""NLFM信号生成"""
xt = self.hamming_nlfm_generate() # 信号：e^(j*phase_w)
bool_t = (np.abs(td[n, :]) < self.PW / 2) #当前时间在脉冲持续时间内
coords = np.where(bool_t)
start_index = coords[0][0]
end_index = start_index + xt.shape[0]
baseband_signal = bool_t.astype(complex)
baseband_signal[start_index:end_index] = baseband_signal[start_index:end_index] * xt
echo_signal[n, :] = (self.RCS * self.PA * baseband_signal * np.exp(
1j * 2 * np.pi * self.f0 * td[n, :] + 1j * self.theta0))
elif signal_type == 5:
xt = self.costas_signal_generate()
bool_t = (np.abs(td[n, :]) < self.PW / 2)
coords = np.where(bool_t)
start_index = coords[0][0]
end_index = start_index + xt.shape[0]
baseband_signal = bool_t.astype(complex)
baseband_signal[start_index:end_index] = baseband_signal[start_index:end_index] * xt
echo_signal[n, :] = (self.RCS * self.PA * baseband_signal * np.exp(
1j * 2 * np.pi * self.f0 * td[n, :] + 1j * self.theta0))
""" 时域波形 """
# plt.figure(figsize=(8, 6))
# plt.plot(tr[0, :], echo_signal[0, :])
# # plt.title("加入自由空间路损后的快时间信号")
# plt.xlabel("时间/s")
# plt.ylabel("幅度")
# plt.grid()
# plt.show()
""" 频域分析 """
# plt.figure(figsize=(10, 8))
# plt.plot(np.abs(np.fft.fft(echo_signal[0, :])))
# plt.grid()
# plt.show()
if signal_type == 1 or signal_type == 4:
return tr, td, echo_signal
elif signal_type == 2 or signal_type == 3 or signal_type == 5:
return tr, td, echo_signal
def add_noise(self, signal):
"""
添加高斯白噪声
"""
noise_signal = np.zeros_like(signal)
for n in range(signal.shape[0]):
P_signal = np.mean(np.abs(signal[n, :]) ** 2)
noise_power = P_signal / 10 ** (self.SNR / 10)
noise_signal[n, :] = np.random.normal(0, np.sqrt(noise_power), signal[n, :].shape)
# P_noise = np.mean(np.abs(noise_signal) ** 2)
# P_signal = np.mean(np.abs(signal) ** 2)
# SNR_exp = 10*np.log10(P_signal / P_noise)
return signal + noise_signal
def add_pathlose(self, signal, path, v):
"""
添加自由空间衰减
"""
signal_pathlose = np.zeros_like(signal)
for n in range(signal.shape[0]):
alpha = self.lambdas / (4 * np.pi * 2 * (path - v * n * self.pri))
signal_pathlose[n, :] = alpha * signal[n, :]
return signal_pathlose
def down_converision(self, t, signal, lo_freq):
"""
数字下变频
"""
# if signal_type == 1 or signal_type == 2 or signal_type == 4:
signal_down = signal * np.exp(-1j * 2 * np.pi * lo_freq * t)
# elif signal_type == 3:
# """ """
return signal_down
def code_sampling(bits, N, bw_of_signal, fs):
"""
输入基带信号生成
"""
original_list = bits
# 目标长度
n_target = N
# 码元宽度,基带信号码元宽度等于带宽的倒数
symbol_t = 1/bw_of_signal
# 一个码元对应的采样点数
symbol_N = int(np.around(symbol_t * fs))
# 目标长度对应的码元数量
symbol_num = int(n_target / symbol_N)
# 完整序列数量以及不完整序列位数
seq_all_num = symbol_num // len(original_list)
seq_residual = symbol_num % len(original_list)
# 每个元素的采样数量
samples_per_element = n_target // len(original_list)
# 扩展后的列表
expanded_list = []
for i in range(seq_all_num):
for element in original_list:
expanded_list.extend([element] * symbol_N)
for element in original_list[0:seq_residual]:
expanded_list.extend([element] * symbol_N)
# 如果列表长度小于目标长度，补齐最后一个元素
while len(expanded_list) < n_target:
expanded_list.append(original_list[-1])
return expanded_list
def pulse_compression(signal, reference_signal, PRT_num, Nwindow, N):
"""
PC脉冲压缩: 对接收回波进行匹配滤波
Param:
signal: 输入的回波信号阵列
reference_signal: 参考信号
PRT_num: 脉冲数量
Nwindow: 回波窗内采样点数
N: 参考信号采样点数
"""
Nfft = 2 ** int(np.ceil(np.log2(Nwindow + N - 1))) # 方便使用FFT算法，满足2的次方形式
H_filter_w = np.zeros((PRT_num, Nfft), dtype=complex)
Srw = np.zeros((PRT_num, Nfft), dtype=complex)
for n in range(PRT_num):
H_filter_w[n, :] = np.fft.fftshift(np.fft.fft(reference_signal, Nfft))
Srw[n, :] = np.fft.fftshift(np.fft.fft(signal[n, :], Nfft))
S_output_w = H_filter_w * Srw
s_output_ifft = np.fft.ifft(S_output_w, Nfft, axis=1) # 时域压缩后的结果（按行进行ifft）
s_output_t = s_output_ifft[:, :N + Nwindow - 1] # 每一行取 (1:N+Nwind-1) 个点即为线性卷积的结果
return s_output_t
def MTI_processing(signal, t):
"""
MTI处理：单延迟线对消器时域实现
Param:
signal: 输入的接收信号阵列 (二维数组: 速度-距离)
t: 回波窗时间轴(快时间时间轴)
"""
# 创建存储MTI后结果的矩阵,横轴为慢时间，纵轴为快时间
mti_signal = np.zeros([signal.shape[0] - 1, signal.shape[1]], dtype=complex)
for i in range(1, signal.shape[0] - 1):
mti_signal[i-1, :] = signal[i, :] - signal[i-1, :]
""" 频域分析 """
# plt.figure(figsize=(8, 6))
# plt.plot(t, np.abs(mti_signal[0, :]))
# # plt.plot(t, np.abs(signal[1, :]))
# plt.grid()
# plt.tight_layout()
# # plt.show()
# 距离轴
distance = t * C / 2
# 脉冲编号
pulse_list = np.linspace(0, mti_signal.shape[0], mti_signal.shape[0])
plt.figure(figsize=(8, 6))
plt.imshow(np.abs(mti_signal), cmap='jet', interpolation='nearest',
aspect='auto', extent=[distance[0], distance[-1],
pulse_list[-1], pulse_list[0]])
plt.ylabel('脉冲ID号')
plt.xlabel('距离(m)')
plt.title('脉冲-距离 MTI结果图')
# plt.grid()
plt.tight_layout()
# plt.show()
return mti_signal
def MTD_processing(mti_signal, t, pri, fc, NFFT):
"""
MTD处理：实现快慢时间2维FFT，在慢时间维度进行FFT
Param:
mti_signal: 输入的MTI信号阵列 (二维数组: 快拍-距离)
t: 回波窗时间轴
pri: 脉冲重复间隔
fc: 载频频率
"""
# 进行慢时间FFT
prf = 1 / pri
# Nfft_slowtime = int(prf)
Nfft_slowtime = NFFT
slow_time_fft = np.fft.fftshift(np.fft.fft(mti_signal, Nfft_slowtime, axis=0), axes=0)
# 获取距离和速度轴
distance = t * C / 2
# 绘制距离维FFT图谱
# 脉冲编号
PRF_space = np.linspace(-prf / 2, prf / 2, Nfft_slowtime)
V_space = PRF_space * (C / fc) / 2
plt.figure(figsize=(8, 6))
plt.imshow(np.abs(slow_time_fft), cmap='jet', interpolation='nearest',
aspect='auto', extent=[distance[0], distance[-1], V_space[-1], V_space[0]])
plt.colorbar(label='Amplitude')
plt.ylabel('速度(m/s)')
plt.xlabel('距离(m)')
plt.title('速度-距离 MTD结果图')
plt.tight_layout()
return np.abs(slow_time_fft), V_space
def CFAR_processing(mtd_signal, CFAR_type, num_protectUnit, num_referUnit, CFAR_gain, N_Range):
# """
# 恒虚警处理: 信号检测门限
# """
num_protect_unit = num_protectUnit * N_Range
num_refer_unit = num_referUnit * N_Range
len_signal = mtd_signal.shape[1]
threshold_cfar = np.zeros_like(mtd_signal)
detection_result = np.zeros_like(mtd_signal)
target_num = np.zeros((mtd_signal.shape[0], 1), dtype=int)
for n in range(mtd_signal.shape[0]):
for index_signal in range(len_signal):
average_l = 0 # Mean of the left reference window
average_r = 0 # Mean of the right reference window
ref_l1 = index_signal - num_protect_unit - num_refer_unit # Left reference window start左参考单元起点
ref_l2 = index_signal - num_protect_unit - 1 # Left reference window end左参考单元终点
ref_r1 = index_signal + num_protect_unit + 1 # Right reference window start
ref_r2 = index_signal + num_protect_unit + num_refer_unit # Right reference window end
# Calculate mean energy of reference windows
if ref_l1 >= 0:
average_l = np.mean(mtd_signal[n, ref_l1:ref_l2 + 1])
if ref_r2 < len_signal:
average_r = np.mean(mtd_signal[n, ref_r1:ref_r2 + 1])
if CFAR_type == 1:
# GO_CFAR
ref_average_max = max(average_l, average_r)
elif CFAR_type == 2:
# SO_CFAR
ref_average_max = min(average_l, average_r)
elif CFAR_type == 3:
# CA_CFAR
ref_average_max = (average_l + average_r) / 2
else:
raise ValueError("Invalid CFAR type. Use 1 for GO_CFAR, 2 for SO_CFAR, or 3 for CA_CFAR.")
threshold_cfar[n, index_signal] = ref_average_max * CFAR_gain # Calculate CFAR threshold
""" 先判断该位置处的信号是否超过检测门限，若超过检测门限则判断其是否为峰值点 """
if mtd_signal[n, index_signal] > threshold_cfar[n, index_signal]:
if 0 < index_signal < len_signal - 1:
if (mtd_signal[n, index_signal] > mtd_signal[n, index_signal - 1] and mtd_signal[n, index_signal] >
mtd_signal[n, index_signal + 1]):
# Peak point detection
target_num[n] += 1
detection_result[n, index_signal] = mtd_signal[n, index_signal]
return detection_result, threshold_cfar, target_num
帮我解释函数的输入输出以及实现功能
最新发布
扫一扫
3744
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
