791 lines
32 KiB
Plaintext
791 lines
32 KiB
Plaintext
# HRRP（高分辨率距离像算法）
|
||
|
||
## 算法概述
|
||
|
||
高分辨率距离像（High Resolution Range Profile, HRRP）是雷达信号处理中的核心技术，通过对雷达回波信号进行脉冲压缩和多普勒处理，获取目标在距离维上的散射点分布特征。HRRP具有以下特点：
|
||
|
||
- 距离高分辨：能够区分目标上相邻较近的散射中心
|
||
- 特征稳定：为目标识别提供物理特征描述
|
||
- 计算高效：相比二维成像计算复杂度较低
|
||
- 实时性强：适用于实时目标识别系统
|
||
|
||
## 代码来源
|
||
<!-- - [Keystone.cpp](../../_static/hrrp/Keystone.cpp)
|
||
- [FFT4f.cpp](../../_static/hrrp/FFT4f.cpp)
|
||
- [FFT4f.h](../../_static/hrrp/FFT4f.h) -->
|
||
|
||
## MindSpore Signal+ 实现
|
||
将原有代码重构为基于计算图的MindSpore Cell，实现雷达信号处理流程模块化。
|
||
|
||
### 1 数据读取
|
||
在Python中使用MindSpore Signal+时，我们可以使用NumPy的fromfile接口读二进制格式文件，因此在Python代码开头需要导入NumPy库。
|
||
|
||
示例代码：
|
||
```python
|
||
def read_hrrp_data(filename, maichongshu, jln, rfftn):
|
||
"""
|
||
读取HRRP二进制数据文件并处理
|
||
"""
|
||
with open(filename, 'rb') as fp:
|
||
expected_size = maichongshu * jln * 2
|
||
data = np.fromfile(fp, dtype=np.float32)
|
||
if len(data) != expected_size:
|
||
raise ValueError(f"数据长度错误：期望{expected_size}个点，实际读取{len(data)}个点")
|
||
data = data.reshape(maichongshu, jln, 2)
|
||
datar = data[:, :, 0]
|
||
datai = data[:, :, 1]
|
||
# 创建填充0的完整数组
|
||
datar_full = np.zeros((maichongshu, rfftn), dtype=np.float32)
|
||
datai_full = np.zeros((maichongshu, rfftn), dtype=np.float32)
|
||
datar_full[:, :jln] = datar
|
||
datai_full[:, :jln] = datai
|
||
return datar_full, datai_full
|
||
```
|
||
|
||
读到的数据需要转换为复数形式，以便后续处理。
|
||
```python
|
||
datar, datai = read_hrrp_data("hrrp_400_512.dat", maichongshu, jln, rfftn)
|
||
input_data = np.zeros(datar.shape, dtype=np.complex64)
|
||
input_data = datar + 1j * datai
|
||
```
|
||
|
||
### 2 数据预处理
|
||
从算法整体分析，数据读取后到核心计算之前的步骤，主要是对相位补偿因子计算，Keystone变换参数预计算，泰勒窗函数预计算，循环移位索引预计算。这部分都是核心计算的前期准备，建议将这部分代码封装在`__init__`函数中完成，不放在`construct`函数中可以避免额外的开销，当实例化一个类时自动触发一次`__init__`函数。
|
||
示例代码：
|
||
```python
|
||
class Hrrp(nn.Cell):
|
||
def __init__(self, hangmc, liejl, f0, fs, B, mohuhalfshu, PRF, vdengxiao, mubiaojuli):
|
||
super(Hrrp, self).__init__()
|
||
self.hangmc = hangmc
|
||
self.liejl = liejl
|
||
self.f0 = f0
|
||
self.fs = fs
|
||
self.B = B
|
||
self.mohuhalfshu = mohuhalfshu
|
||
self.PRF = PRF
|
||
self.vdengxiao = vdengxiao
|
||
self.mubiaojuli = mubiaojuli
|
||
self.fft = mr.FFT()
|
||
self.ifft = mr.IFFT()
|
||
self.fftshift = mr.FFTShift()
|
||
self.abs = mr.ComplexAbs()
|
||
self.stoltsun = mr.Stoltsun(dim=0)
|
||
self.mul = ops.Mul()
|
||
self.exp = ops.Exp()
|
||
self.reduce_sum = ops.ReduceSum(keep_dims=True)
|
||
self.gather = ops.Gather()
|
||
wl = LC / self.f0
|
||
Ka = 2.0 * self.vdengxiao**2 / (wl * self.mubiaojuli)
|
||
mctime = (np.arange(self.hangmc) - self.hangmc * 0.5) / self.PRF
|
||
xiangwei = np.pi * Ka * mctime**2
|
||
self.phase_comp = np.exp(1j * xiangwei).astype(np.complex64)
|
||
self.phase_comp = self.phase_comp.reshape(self.hangmc, 1)
|
||
self.phase_comp = ms.Tensor(self.phase_comp)
|
||
fwn = self.hangmc
|
||
self.fwn = fwn
|
||
rnfft = Get2intm(self.liejl)
|
||
self.rnfft = rnfft
|
||
fwn_8 = (fwn + 7) // 8 * 8 # 8字节对齐
|
||
self.fwn_8 = fwn_8
|
||
self.xold = (np.arange(fwn_8, dtype=np.float32) - fwn / 2.0)[:, np.newaxis]
|
||
self.xold = np.broadcast_to(self.xold, (fwn_8, rnfft))
|
||
i_arr = np.arange(rnfft, dtype=np.float32)
|
||
sigma_arr = (f0 + (i_arr - rnfft * 0.5) * fs / rnfft) / f0
|
||
self.xnew = self.xold * sigma_arr
|
||
self.xold = ms.Tensor(self.xold)
|
||
self.xnew = ms.Tensor(self.xnew)
|
||
self.mchtr = self.TaiLeWindow(self.hangmc, 50)
|
||
window_2d = self.mchtr[:, np.newaxis] # 转换为列向量便于广播
|
||
self.window_2d = ms.Tensor(window_2d, dtype=ms.float32)
|
||
self.j_constant = ms.Parameter(ms.Tensor(1.0j, dtype=ms.complex64), name="j_constant")
|
||
self._precompute_cirshift_indices()
|
||
|
||
def _precompute_cirshift_indices(self):
|
||
rnfft = self.rnfft
|
||
shift = self.rnfft // 2
|
||
indices = np.arange(rnfft, dtype=np.int32)
|
||
indices_shifted = (indices - shift) % rnfft
|
||
self.cirshift_indices = ms.Tensor(indices_shifted, dtype=ms.int32)
|
||
|
||
|
||
def construct(self, input_data):
|
||
...
|
||
```
|
||
|
||
#### 2.1 相位补偿因子计算
|
||
```python
|
||
wl = LC / self.f0 # 计算波长
|
||
Ka = 2.0 * self.vdengxiao**2 / (wl * self.mubiaojuli) # 调频斜率
|
||
mctime = (np.arange(self.hangmc) - self.hangmc * 0.5) / self.PRF # 时间轴
|
||
xiangwei = np.pi * Ka * mctime**2 # 相位变化量
|
||
self.phase_comp = np.exp(1j * xiangwei).astype(np.complex64) # 相位补偿因子
|
||
```
|
||
功能说明：
|
||
|
||
- 计算二次相位误差补偿因子，用于校正由于目标运动引起的相位畸变
|
||
- 将连续时间计算转换为离散时间序列，便于数字信号处理
|
||
- 相位补偿因子形状为 (hangmc, 1)，便于广播到所有距离单元
|
||
|
||
#### 2.2 Keystone变换参数预计算
|
||
```python
|
||
# 关键参数计算
|
||
fwn = self.hangmc
|
||
self.fwn = fwn
|
||
rnfft = Get2intm(self.liejl) # 获取2的幂次方长度
|
||
self.rnfft = rnfft
|
||
fwn_8 = (fwn + 7) // 8 * 8 # 8字节对齐优化
|
||
self.fwn_8 = fwn_8
|
||
# 方位向坐标轴生成
|
||
self.xold = (np.arange(fwn_8, dtype=np.float32) - fwn / 2.0)[:, np.newaxis]
|
||
self.xold = np.broadcast_to(self.xold, (fwn_8, rnfft))
|
||
# 距离频率轴生成和缩放因子计算
|
||
i_arr = np.arange(rnfft, dtype=np.float32)
|
||
sigma_arr = (f0 + (i_arr - rnfft * 0.5) * fs / rnfft) / f0
|
||
self.xnew = self.xold * sigma_arr
|
||
# 转换为MindSpore Tensor
|
||
self.xold = ms.Tensor(self.xold)
|
||
self.xnew = ms.Tensor(self.xnew)
|
||
```
|
||
参数设计原理：
|
||
| 参数 | 计算公式 | 物理意义 |
|
||
|-----------| ---------------------- | -------------------- |
|
||
| xold | (i - fwn/2.0) | 方位向原始坐标，零中心化 |
|
||
| sigma_arr | (f0 + (i - rnfft/2)*fs/rnfft)/f0 | 距离徙动校正因子 |
|
||
| xnew | xold * sigma_arr | Keystone变换后的坐标 |
|
||
|
||
|
||
#### 2.3 泰勒窗函数预计算
|
||
```python
|
||
self.mchtr = self.TaiLeWindow(self.hangmc, 50)
|
||
window_2d = self.mchtr[:, np.newaxis]
|
||
self.window_2d = ms.Tensor(window_2d, dtype=ms.float32)
|
||
|
||
def TaiLeWindow(self, N, param2):
|
||
# 参数处理
|
||
N = max(N, 1)
|
||
nLevel = 5
|
||
nsll = 45 if param2 < 13.0 else param2
|
||
# 参数转换
|
||
mB = math.pow(10.0, nsll / 20.0)
|
||
mA = math.log(mB + math.sqrt(mB * mB - 1.0)) / math.pi
|
||
msquaQ = (nLevel * nLevel) / (mA * mA + (nLevel - 0.5) * (nLevel - 0.5))
|
||
# 计算窗函数系数
|
||
m_dFm = [0.0] * nLevel
|
||
for m in range(nLevel - 1):
|
||
mtmpa = 1.0
|
||
mtmpb = 1.0
|
||
for i in range(nLevel - 1):
|
||
# 计算第一个乘积项
|
||
term = 1 - (m+1)**2 / (msquaQ * (mA**2 + (i+0.5)**2))
|
||
mtmpa *= term
|
||
# 跳过相同索引
|
||
if m == i:
|
||
continue
|
||
# 计算第二个乘积项
|
||
term2 = 1 - (m+1)**2 / ((i+1)**2)
|
||
mtmpb *= term2
|
||
# 计算系数值
|
||
m_dFm[m] = (pow(-1, m+2) * mtmpa) / (2.0 * mtmpb)
|
||
# 计算窗函数值
|
||
h = np.zeros(N)
|
||
for i in range(N):
|
||
dwin = 0.0
|
||
for m in range(nLevel - 1):
|
||
angle = 2 * math.pi * (m+1) * (i - N/2.0 + 0.5) / N
|
||
dwin += m_dFm[m] * math.cos(angle)
|
||
h[i] = 1.0 + 2.0 * dwin
|
||
# 归一化处理
|
||
mean_val = np.mean(h)
|
||
if abs(mean_val) > 1e-15:
|
||
h /= mean_val
|
||
return h
|
||
```
|
||
窗函数作用：
|
||
- 减少频谱泄露，提高频率分辨率
|
||
- 抑制旁瓣，提高主瓣分辨率
|
||
- 参数50表示旁瓣抑制比约为50dB
|
||
|
||
#### 2.4 循环移位索引预计算
|
||
```python
|
||
def _precompute_cirshift_indices(self):
|
||
rnfft = self.rnfft
|
||
shift = self.rnfft // 2 # 中心频率点索引
|
||
indices = np.arange(rnfft, dtype=np.int32)
|
||
indices_shifted = (indices - shift) % rnfft
|
||
self.cirshift_indices = ms.Tensor(indices_shifted, dtype=ms.int32)
|
||
```
|
||
循环移位的作用：
|
||
- 频谱中心化：将零频分量移动到频谱中心
|
||
- 视觉优化：便于频谱显示和分析
|
||
- 算法要求：某些算法（如IFFT）要求输入是中心对称的
|
||
|
||
|
||
### 3 核心计算
|
||
|
||
步骤分解：
|
||
1. 相位补偿：首先对输入数据（复数形式）进行相位补偿，乘以预计算好的相位补偿因子（用于校正二次相位误差）。
|
||
2. 多模糊数处理：计算需要搜索的模糊数总数（mohushu = 2 * mohuhalfshu + 1），然后对每一个模糊数进行以下操作：
|
||
|
||
a. 计算当前模糊数对应的偏移量（offset = mh - mohuhalfshu），这个偏移量将用于Keystone变换中的相位补偿。
|
||
|
||
b. 进行Keystone变换（Keystone_interp），这是HRRP算法的核心，用于校正距离徙动。
|
||
|
||
c. 对Keystone变换后的数据，进行HRRP生成（Gethrrp），得到每一距离单元的幅度和相位信息，同时返回处理后的数据（result）和幅度数据（hrrp）。
|
||
|
||
d. 计算该HRRP的熵（Gethrrp_v2shang），熵值越小表示能量越集中，成像质量越好。
|
||
|
||
e. 将当前模糊数对应的结果（result）和熵值（shang_value）分别保存到列表results和shang中。
|
||
|
||
3. 选择最优结果：将所有模糊数对应的结果堆叠起来（使用ops.cat），然后找到熵最小的那个索引（xuhao = ops.argmin(shang)），并取出该索引对应的结果（min_result = results[xuhao]）。
|
||
|
||
4. 更新输出：将最优结果的第一行（即HRRP处理后的第一行，包含了最优的相位和幅度信息）与原始输入数据的第二行到最后一行拼接起来，形成最终的输出。这里之所以只替换第一行，是因为在HRRP处理中，我们只关心每个距离单元的最强散射点，而原始数据的第一行被用来存储这个信息。
|
||
|
||
核心算法整体流程：
|
||
```python
|
||
class Hrrp(nn.Cell):
|
||
def __init__(self, hangmc, liejl, f0, fs, B, mohuhalfshu, PRF, vdengxiao, mubiaojuli):
|
||
super(Hrrp, self).__init__()
|
||
self.hangmc = hangmc
|
||
self.f0 = f0
|
||
self.fs = fs
|
||
self.B = B
|
||
self.mohuhalfshu = mohuhalfshu
|
||
self.PRF = PRF
|
||
self.vdengxiao = vdengxiao
|
||
self.mubiaojuli = mubiaojuli
|
||
...
|
||
wl = LC / self.f0
|
||
Ka = 2.0 * self.vdengxiao**2 / (wl * self.mubiaojuli)
|
||
mctime = (np.arange(self.hangmc) - self.hangmc * 0.5) / self.PRF
|
||
xiangwei = np.pi * Ka * mctime**2
|
||
self.phase_comp = np.exp(1j * xiangwei).astype(np.complex64)
|
||
self.phase_comp = self.phase_comp.reshape(self.hangmc, 1)
|
||
self.phase_comp = ms.Tensor(self.phase_comp)
|
||
...
|
||
|
||
def construct(self, input_data):
|
||
out = input_data * self.phase_comp
|
||
mohushu = 2 * self.mohuhalfshu + 1
|
||
results = []
|
||
shang = []
|
||
for mh in range(mohushu):
|
||
offset = mh - self.mohuhalfshu
|
||
tmp = self.Keystone_interp(out, self.f0, self.fs, self.B, offset)
|
||
hrrp, result = self.Gethrrp(tmp, self.hangmc, self.liejl)
|
||
shang_value = self.Gethrrp_v2shang(hrrp)
|
||
shang.append(shang_value)
|
||
result = ops.expand_dims(result, 0)
|
||
results.append(result)
|
||
results = ops.cat(results, axis=0)
|
||
shang = ops.cat(shang)
|
||
xuhao = ops.argmin(shang)
|
||
min_result = results[xuhao]
|
||
out = ops.cat([min_result[0:1, :], out[1:, :]], axis=0)
|
||
return out
|
||
```
|
||
|
||
#### 3.1 相位补偿
|
||
```python
|
||
out = input_data * self.phase_comp
|
||
```
|
||
物理意义：
|
||
- 补偿目标：校正由于雷达与目标相对运动引起的二次相位误差
|
||
- 误差来源：目标斜距变化引起的时间延迟变化
|
||
- 数学形式：补偿因子为 exp(j·π·Ka·t²)，其中Ka为调频斜率
|
||
|
||
#### 3.2 多模糊数处理循环
|
||
1. 模糊数范围确定：
|
||
```python
|
||
mohushu = 2 * self.mohuhalfshu + 1
|
||
```
|
||
设计原理：
|
||
- 由于多普勒频率模糊，真实的调频斜率有多个候选值
|
||
- 搜索范围：[-mohuhalfshu, +mohuhalfshu]
|
||
- 典型设置：mohuhalfshu=5 → 搜索11个候选值
|
||
|
||
2. 循环处理每个模糊数：
|
||
```python
|
||
for mh in range(mohushu):
|
||
offset = mh - self.mohuhalfshu # 当前模糊数偏移
|
||
tmp = self.Keystone_interp(out, self.f0, self.fs, self.B, offset)
|
||
hrrp, result = self.Gethrrp(tmp, self.hangmc, self.liejl)
|
||
shang_value = self.Gethrrp_v2shang(hrrp)
|
||
shang.append(shang_value)
|
||
result = ops.expand_dims(result, 0)
|
||
results.append(result)
|
||
```
|
||
处理流程分解：
|
||
|步骤 |功能 | 关键操作 |
|
||
|--------------| ----------- | -------------------- |
|
||
|Keystone变换 |距离徙动校正 | 使用当前模糊数进行插值 |
|
||
|HRRP生成 |距离像提取 | 方位向FFT，提取峰值 |
|
||
|熵计算 | 质量评估 | 计算HRRP的熵值 |
|
||
|
||
#### 3.3 Keystone变换详解
|
||
1. 变换原理
|
||
Keystone变换是一种距离徙动校正算法，用于解决以下问题：
|
||
<div align="center">
|
||
<img src="../../_static/Keystone.png" width="50%">
|
||
</div>
|
||
|
||
变换公式：
|
||
原始坐标：t (方位时间), f (距离频率)
|
||
变换后：t' = t × (f₀/(f₀+f))
|
||
|
||
2. 代码实现
|
||
```python
|
||
# 1. 距离向FFT（转到距离频率域）
|
||
# 2. 相位补偿（补偿距离徙动），方位向插值（Keystone变换核心）
|
||
# 3. 距离向IFFT（转回时域）
|
||
def Keystone_interp(self, input_data, f0, fs, B, moshuzhouqishu):
|
||
# 获取维度
|
||
fwn = self.fwn
|
||
fwn_8 = self.fwn_8
|
||
# 计算距离FFT长度
|
||
rnfft = self.rnfft
|
||
# 模块1: 距离向FFT
|
||
out = self.FFT4f(input_data, 1)
|
||
nyquist_idx = rnfft//2
|
||
new_nyquist_values = (out[:, nyquist_idx-1] + out[:, nyquist_idx+1]) / 2.0
|
||
# 使用concat操作替换特定列
|
||
left_part = out[:, :nyquist_idx]
|
||
nyquist_part = new_nyquist_values.expand_dims(1) # 增加维度
|
||
right_part = out[:, nyquist_idx+1:]
|
||
# 重新拼接
|
||
out = ops.concat([left_part, nyquist_part, right_part], axis=1)
|
||
# 计算无小点数量
|
||
wuxiaodiansujl = int(np.floor(0.5 * (fs - B) * rnfft / fs))
|
||
wuxiaodiansujl = (wuxiaodiansujl >> 3) << 3 # 8字节对齐
|
||
wuxiaodiansujl = max(0, wuxiaodiansujl)
|
||
# 模块2: 方位向插值
|
||
pindian = (ops.arange(0, rnfft, dtype=ms.float32) - 0.5 * rnfft) * fs / rnfft
|
||
j_indices = ops.arange(0, fwn, dtype=ms.float32).reshape(-1, 1)
|
||
angle = -2 * PI * j_indices * pindian / f0 * moshuzhouqishu
|
||
complex_angle = angle.astype(ms.complex64)
|
||
j_times_angle = complex_angle * self.j_constant
|
||
phase_comp_all = self.exp(j_times_angle)
|
||
temp_complex_top = out * phase_comp_all
|
||
temp_complex = ops.Pad(((0, fwn_8 - fwn), (0, 0)))(temp_complex_top)
|
||
temp_complex_valid = temp_complex[:, wuxiaodiansujl:rnfft - wuxiaodiansujl]
|
||
xnew_valid = self.xnew[:, wuxiaodiansujl:rnfft - wuxiaodiansujl]
|
||
xold_valid = self.xold[:, wuxiaodiansujl:rnfft - wuxiaodiansujl]
|
||
tmp = self.stoltsun(temp_complex_valid, xnew_valid, xold_valid)
|
||
update_part = tmp
|
||
left_part = out[:fwn, :wuxiaodiansujl]
|
||
right_part = out[:fwn, rnfft - wuxiaodiansujl:]
|
||
# 在列方向拼接三部分
|
||
out = ops.concat([left_part, update_part, right_part], axis=1)
|
||
# 模块3: 距离向IFFT
|
||
out = self.cirshift_optimized(out)
|
||
out = self.IFFT4f(out, 0)
|
||
return out
|
||
```
|
||
|
||
#### 3.4 HRRP生成
|
||
```python
|
||
# 1. 方位向加窗（泰勒窗）
|
||
# 2. 方位向FFT（转到多普勒域）
|
||
# 3. 峰值提取（每个距离单元的幅度）
|
||
# 4. 返回HRRP幅度和复数值
|
||
def Gethrrp(self, input_data, hangmc, liejln):
|
||
afftn = Get2intm(hangmc)
|
||
# 只取有效数据部分进行加窗处理
|
||
windowed_data = input_data[:hangmc, :liejln] * self.window_2d
|
||
windowed_data = windowed_data.astype(ms.complex64)
|
||
windowed_data_t = ops.transpose(windowed_data, (1, 0))
|
||
padding = ((0, 0), (0, afftn - windowed_data_t.shape[1]))
|
||
fft_input = ops.Pad(paddings=padding)(windowed_data_t)
|
||
fft_ret = self.FFT4f(fft_input, 0)
|
||
# 计算幅度平方 - 向量化操作
|
||
# fd = self.abs(fft_ret)**2
|
||
abs_result = ops.cast(ops.abs(fft_ret), ms.float32)
|
||
fd = abs_result**2
|
||
# 找每列的峰值位置 - 向量化操作
|
||
maxp_indices = ops.argmax(fd, dim=1)
|
||
# 提取峰值幅度 - 向量化操作
|
||
row_indices = ops.arange(liejln, dtype=ms.int32)
|
||
hrrp = fd[row_indices, maxp_indices]
|
||
peak_values = fft_ret[row_indices, maxp_indices]
|
||
peak_values_2d = ops.expand_dims(peak_values, 0)
|
||
result = ops.cat([peak_values_2d, input_data[1:, :]], axis=0)
|
||
return hrrp, result
|
||
```
|
||
|
||
HRRP定义：
|
||
- 高分辨率距离像 = 目标在距离维的散射点分布
|
||
- 物理意义：目标的"电磁指纹"，用于目标识别
|
||
|
||
#### 3.5 熵值评估
|
||
```python
|
||
def Gethrrp_v2shang(self, x):
|
||
sum_sq = self.reduce_sum(x**2) # 总能量
|
||
is_nan = sum_sq != sum_sq
|
||
# 检查无效情况
|
||
is_valid = ops.logical_and(sum_sq > 1e-10, ops.logical_not(is_nan))
|
||
# 归一化并计算熵
|
||
t = x / (sum_sq)
|
||
log_t = ops.log(t)
|
||
entropy = -self.reduce_sum(t * log_t) # 香农熵
|
||
entropy = entropy.reshape(1)
|
||
# 如果无效则返回0
|
||
result = ops.where(is_valid, entropy, ops.zeros(1, ms.float32))
|
||
result = result.reshape(1)
|
||
return result
|
||
```
|
||
|
||
熵的物理意义：
|
||
|HRRP特征 |熵值表现 |物理解释|
|
||
|----------|------------|--------|
|
||
|能量集中 |低熵值 |少数强散射点，成像清晰|
|
||
|能量分散 |高熵值 |多个弱散射点，成像模糊|
|
||
|噪声影响 |高熵值 |背景噪声使能量分散|
|
||
|
||
熵值选择准则：熵值越小 → 能量越集中 → 成像质量越好 → 模糊数越准确
|
||
|
||
#### 3.6 最优结果选择
|
||
```python
|
||
# 合并所有结果
|
||
results = ops.cat(results, axis=0) # [mohushu, hangmc, liejl]
|
||
shang = ops.cat(shang) # [mohushu]
|
||
# 找到最小熵对应的索引
|
||
xuhao = ops.argmin(shang) # 最优模糊数索引
|
||
# 选择最优结果
|
||
min_result = results[xuhao] # 最优HRRP结果
|
||
out = ops.cat([min_result[0:1, :], out[1:, :]], axis=0) # 更新最终结果
|
||
```
|
||
|
||
### 4 数据后处理
|
||
这里将结果导出成二进制数据的dat文件，用于和正确结果比较。确认结果无误后，可通过 ``mindspore.export`` 导出 MINDIR 模型，便于在 MindSpore Lite 端部署（板卡侧运行）。
|
||
|
||
## 板卡部署
|
||
模型部署建议使用 ``YHFT-IDE``，它集成了模型转换、模型可视化与 MindSpore Lite 端部署模板。具体使用方法可参考 {ref}`HelloDSP MindSpore Lite端 <c++>`。
|
||
|
||
以下是完整的MindSpore Signal+实现的代码：
|
||
```python
|
||
import mindspore as ms
|
||
import numpy as np
|
||
from mindspore import nn, ops
|
||
import math
|
||
import mindradar as mr
|
||
|
||
def read_hrrp_data(filename, maichongshu, jln, rfftn):
|
||
# 打开二进制文件并读取所有数据
|
||
with open(filename, 'rb') as fp:
|
||
# 计算预期数据量并读取
|
||
expected_size = maichongshu * jln * 2 # 每个复数点包含实部和虚部
|
||
data = np.fromfile(fp, dtype=np.float32)
|
||
# 验证数据完整性
|
||
if len(data) != expected_size:
|
||
raise ValueError(f"数据长度错误：期望{expected_size}个点，实际读取{len(data)}个点")
|
||
|
||
# 重塑数据为3D数组 [脉冲索引, 数据点, 实部/虚部]
|
||
data = data.reshape(maichongshu, jln, 2)
|
||
# 分离实部和虚部
|
||
datar = data[:, :, 0] # 实部 [maichongshu, jln]
|
||
datai = data[:, :, 1] # 虚部 [maichongshu, jln]
|
||
# 创建填充0的完整数组
|
||
datar_full = np.zeros((maichongshu, rfftn), dtype=np.float32)
|
||
datai_full = np.zeros((maichongshu, rfftn), dtype=np.float32)
|
||
# 复制有效数据部分
|
||
datar_full[:, :jln] = datar
|
||
datai_full[:, :jln] = datai
|
||
return datar_full, datai_full
|
||
|
||
def Get2intm(intnumb):
|
||
# 验证输入有效性
|
||
if not isinstance(intnumb, int) or intnumb <= 0:
|
||
raise ValueError("Input must be a positive integer")
|
||
# 特殊处理1的情况（log2(1)=0）
|
||
if intnumb == 1:
|
||
return 1
|
||
# 使用更精确的log2计算
|
||
log2_val = math.log2(intnumb)
|
||
# 向上取整获取指数
|
||
nn = math.ceil(log2_val)
|
||
# 计算2的nn次方
|
||
result = 2 ** nn
|
||
return int(result) # 转换为整数类型
|
||
|
||
|
||
def write_hrrp_binary_v2(filename, datar, datai, maichongshu, jln):
|
||
# 参数强验证
|
||
if datar.shape != (maichongshu, 512) or datai.shape != (maichongshu, 512):
|
||
raise ValueError(f"数组维度必须为 ({maichongshu}, 512)")
|
||
# 创建交错存储矩阵
|
||
interleaved = np.empty((maichongshu, jln*2), dtype=np.float32)
|
||
# 精确截取前2560列
|
||
interleaved[:, ::2] = datar[:, :jln] # 实部
|
||
interleaved[:, 1::2] = datai[:, :jln] # 虚部
|
||
# 原子写入操作
|
||
interleaved.tofile(filename)
|
||
|
||
class Hrrp(nn.Cell):
|
||
def __init__(self, hangmc, liejl, f0, fs, B, mohuhalfshu, PRF, vdengxiao, mubiaojuli):
|
||
super(Hrrp, self).__init__()
|
||
self.hangmc = hangmc
|
||
self.liejl = liejl
|
||
self.f0 = f0
|
||
self.fs = fs
|
||
self.B = B
|
||
self.mohuhalfshu = mohuhalfshu
|
||
self.PRF = PRF
|
||
self.vdengxiao = vdengxiao
|
||
self.mubiaojuli = mubiaojuli
|
||
self.fft = mr.FFT()
|
||
self.ifft = mr.IFFT()
|
||
self.fftshift = mr.FFTShift()
|
||
self.stoltsun = mr.Stoltsun(dim=0)
|
||
self.mul = ops.Mul()
|
||
self.exp = ops.Exp()
|
||
self.reduce_sum = ops.ReduceSum(keep_dims=True)
|
||
self.gather = ops.Gather()
|
||
wl = LC / self.f0
|
||
Ka = 2.0 * self.vdengxiao**2 / (wl * self.mubiaojuli)
|
||
mctime = (np.arange(self.hangmc) - self.hangmc * 0.5) / self.PRF
|
||
xiangwei = np.pi * Ka * mctime**2
|
||
self.phase_comp = np.exp(1j * xiangwei).astype(np.complex64)
|
||
self.phase_comp = self.phase_comp.reshape(self.hangmc, 1)
|
||
self.phase_comp = ms.Tensor(self.phase_comp)
|
||
fwn = self.hangmc
|
||
self.fwn = fwn
|
||
rnfft = Get2intm(self.liejl)
|
||
self.rnfft = rnfft
|
||
fwn_8 = (fwn + 7) // 8 * 8 # 8字节对齐
|
||
self.fwn_8 = fwn_8
|
||
self.xold = (np.arange(fwn_8, dtype=np.float32) - fwn / 2.0)[:, np.newaxis]
|
||
self.xold = np.broadcast_to(self.xold, (fwn_8, rnfft))
|
||
i_arr = np.arange(rnfft, dtype=np.float32)
|
||
sigma_arr = (f0 + (i_arr - rnfft * 0.5) * fs / rnfft) / f0
|
||
self.xnew = self.xold * sigma_arr
|
||
self.xold = ms.Tensor(self.xold)
|
||
self.xnew = ms.Tensor(self.xnew)
|
||
self.mchtr = self.TaiLeWindow(self.hangmc, 50)
|
||
window_2d = self.mchtr[:, np.newaxis] # 转换为列向量便于广播
|
||
self.window_2d = ms.Tensor(window_2d, dtype=ms.float32)
|
||
self.j_constant = ms.Parameter(ms.Tensor(1.0j, dtype=ms.complex64), name="j_constant")
|
||
self._precompute_cirshift_indices()
|
||
|
||
def _precompute_cirshift_indices(self):
|
||
rnfft = self.rnfft
|
||
shift = self.rnfft // 2
|
||
# 创建索引数组 [0, 1, 2, ..., rnfft-1]
|
||
indices = np.arange(rnfft, dtype=np.int32)
|
||
# 计算移位后的索引
|
||
# 对于左移shift：new_index = (old_index - shift) mod rnfft
|
||
indices_shifted = (indices - shift) % rnfft
|
||
# 转换为Tensor
|
||
self.cirshift_indices = ms.Tensor(indices_shifted, dtype=ms.int32)
|
||
|
||
|
||
def construct(self, input_data):
|
||
out = input_data * self.phase_comp
|
||
mohushu = 2 * self.mohuhalfshu + 1
|
||
results = []
|
||
shang = []
|
||
for mh in range(mohushu):
|
||
offset = mh - self.mohuhalfshu
|
||
tmp = self.Keystone_interp(out, self.f0, self.fs, self.B, offset)
|
||
hrrp, result = self.Gethrrp(tmp, self.hangmc, self.liejl)
|
||
shang_value = self.Gethrrp_v2shang(hrrp)
|
||
shang.append(shang_value)
|
||
result = ops.expand_dims(result, 0)
|
||
results.append(result)
|
||
results = ops.cat(results, axis=0)
|
||
shang = ops.cat(shang)
|
||
xuhao = ops.argmin(shang)
|
||
min_result = results[xuhao]
|
||
# 将min_result的第一行和out的第二行到最后一行连接起来
|
||
out = ops.cat([min_result[0:1, :], out[1:, :]], axis=0)
|
||
return out
|
||
|
||
def TaiLeWindow(self, N, param2):
|
||
"""
|
||
泰勒窗函数实现
|
||
参数：
|
||
N: 窗口长度
|
||
param2: 加权分贝值 (20~50)
|
||
返回：
|
||
归一化的窗函数数组
|
||
"""
|
||
# 参数处理
|
||
N = max(N, 1)
|
||
nLevel = 5
|
||
nsll = 45 if param2 < 13.0 else param2
|
||
# 参数转换
|
||
mB = math.pow(10.0, nsll / 20.0)
|
||
mA = math.log(mB + math.sqrt(mB * mB - 1.0)) / math.pi
|
||
msquaQ = (nLevel * nLevel) / (mA * mA + (nLevel - 0.5) * (nLevel - 0.5))
|
||
# 计算窗函数系数
|
||
m_dFm = [0.0] * nLevel
|
||
for m in range(nLevel - 1):
|
||
mtmpa = 1.0
|
||
mtmpb = 1.0
|
||
for i in range(nLevel - 1):
|
||
# 计算第一个乘积项
|
||
term = 1 - (m+1)**2 / (msquaQ * (mA**2 + (i+0.5)**2))
|
||
mtmpa *= term
|
||
# 跳过相同索引
|
||
if m == i:
|
||
continue
|
||
# 计算第二个乘积项
|
||
term2 = 1 - (m+1)**2 / ((i+1)**2)
|
||
mtmpb *= term2
|
||
# 计算系数值
|
||
m_dFm[m] = (pow(-1, m+2) * mtmpa) / (2.0 * mtmpb)
|
||
# 计算窗函数值
|
||
h = np.zeros(N)
|
||
for i in range(N):
|
||
dwin = 0.0
|
||
for m in range(nLevel - 1):
|
||
angle = 2 * math.pi * (m+1) * (i - N/2.0 + 0.5) / N
|
||
dwin += m_dFm[m] * math.cos(angle)
|
||
h[i] = 1.0 + 2.0 * dwin
|
||
# 归一化处理
|
||
mean_val = np.mean(h)
|
||
if abs(mean_val) > 1e-15:
|
||
h /= mean_val
|
||
return h
|
||
|
||
def Gethrrp_v2shang(self, x):
|
||
sum_sq = self.reduce_sum(x**2)
|
||
is_nan = sum_sq != sum_sq
|
||
# 检查无效情况
|
||
is_valid = ops.logical_and(sum_sq > 1e-10, ops.logical_not(is_nan))
|
||
# 归一化并计算熵
|
||
t = x / (sum_sq)
|
||
log_t = ops.log(t)
|
||
entropy = -self.reduce_sum(t * log_t)
|
||
entropy = entropy.reshape(1)
|
||
# 如果无效则返回0
|
||
result = ops.where(is_valid, entropy, ops.zeros(1, ms.float32))
|
||
result = result.reshape(1)
|
||
return result
|
||
|
||
def Gethrrp(self, input_data, hangmc, liejln):
|
||
afftn = Get2intm(hangmc)
|
||
# 只取有效数据部分进行加窗处理
|
||
windowed_data = input_data[:hangmc, :liejln] * self.window_2d
|
||
windowed_data = windowed_data.astype(ms.complex64)
|
||
windowed_data_t = ops.transpose(windowed_data, (1, 0))
|
||
padding = ((0, 0), (0, afftn - windowed_data_t.shape[1]))
|
||
fft_input = ops.Pad(paddings=padding)(windowed_data_t)
|
||
fft_ret = self.FFT4f(fft_input, 0)
|
||
# 计算幅度平方 - 向量化操作
|
||
abs_result = ops.cast(ops.abs(fft_ret), ms.float32)
|
||
fd = abs_result**2
|
||
# 找每列的峰值位置 - 向量化操作
|
||
maxp_indices = ops.argmax(fd, dim=1)
|
||
# 提取峰值幅度 - 向量化操作
|
||
row_indices = ops.arange(liejln, dtype=ms.int32)
|
||
hrrp = fd[row_indices, maxp_indices]
|
||
peak_values = fft_ret[row_indices, maxp_indices]
|
||
peak_values_2d = ops.expand_dims(peak_values, 0)
|
||
result = ops.cat([peak_values_2d, input_data[1:, :]], axis=0)
|
||
return hrrp, result
|
||
|
||
def Keystone_interp(self, input_data, f0, fs, B, moshuzhouqishu):
|
||
# 获取维度
|
||
fwn = self.fwn
|
||
fwn_8 = self.fwn_8
|
||
# 计算距离FFT长度
|
||
rnfft = self.rnfft
|
||
# 模块1: 距离向FFT
|
||
out = self.FFT4f(input_data, 1)
|
||
nyquist_idx = rnfft//2
|
||
new_nyquist_values = (out[:, nyquist_idx-1] + out[:, nyquist_idx+1]) / 2.0
|
||
# 使用concat操作替换特定列
|
||
left_part = out[:, :nyquist_idx]
|
||
nyquist_part = new_nyquist_values.expand_dims(1) # 增加维度
|
||
right_part = out[:, nyquist_idx+1:]
|
||
# 重新拼接
|
||
out = ops.concat([left_part, nyquist_part, right_part], axis=1)
|
||
# 计算无小点数量
|
||
wuxiaodiansujl = int(np.floor(0.5 * (fs - B) * rnfft / fs))
|
||
wuxiaodiansujl = (wuxiaodiansujl >> 3) << 3 # 8字节对齐
|
||
wuxiaodiansujl = max(0, wuxiaodiansujl)
|
||
# 模块2: 方位向插值
|
||
pindian = (ops.arange(0, rnfft, dtype=ms.float32) - 0.5 * rnfft) * fs / rnfft
|
||
j_indices = ops.arange(0, fwn, dtype=ms.float32).reshape(-1, 1)
|
||
angle = -2 * PI * j_indices * pindian / f0 * moshuzhouqishu
|
||
complex_angle = angle.astype(ms.complex64)
|
||
j_times_angle = complex_angle * self.j_constant
|
||
phase_comp_all = self.exp(j_times_angle)
|
||
temp_complex_top = out * phase_comp_all
|
||
temp_complex = ops.Pad(((0, fwn_8 - fwn), (0, 0)))(temp_complex_top)
|
||
temp_complex_valid = temp_complex[:, wuxiaodiansujl:rnfft - wuxiaodiansujl]
|
||
xnew_valid = self.xnew[:, wuxiaodiansujl:rnfft - wuxiaodiansujl]
|
||
xold_valid = self.xold[:, wuxiaodiansujl:rnfft - wuxiaodiansujl]
|
||
tmp = self.stoltsun(temp_complex_valid, xnew_valid, xold_valid)
|
||
update_part = tmp
|
||
left_part = out[:fwn, :wuxiaodiansujl]
|
||
right_part = out[:fwn, rnfft - wuxiaodiansujl:]
|
||
# 在列方向拼接三部分
|
||
out = ops.concat([left_part, update_part, right_part], axis=1)
|
||
# 模块3: 距离向IFFT
|
||
# out = self.cirshift(out, rnfft, rnfft // 2)
|
||
out = self.cirshift_optimized(out)
|
||
out = self.IFFT4f(out, 0)
|
||
return out
|
||
|
||
def cirshift_optimized(self, input_data):
|
||
"""优化后的循环移位"""
|
||
# 使用预计算的索引进行gather操作
|
||
return self.gather(input_data, self.cirshift_indices, 1)
|
||
|
||
def cirshift(self, input_data, n, shift):
|
||
shift = shift % n
|
||
if shift == 0:
|
||
return input_data
|
||
if input_data.ndim == 1:
|
||
# 一维情况
|
||
out = ops.concat([input_data[-shift:], input_data[:-shift]])
|
||
else:
|
||
# 二维情况，沿axis=1循环移位
|
||
out = ops.concat([input_data[:, -shift:], input_data[:, :-shift]], axis=1)
|
||
return out
|
||
|
||
def FFT4f(self, input_data, is_inverse):
|
||
fft_ret = self.fft(input_data)
|
||
if is_inverse:
|
||
fft_ret = self.fftshift(fft_ret)
|
||
end_time = time.time() * 1000
|
||
return fft_ret
|
||
|
||
def IFFT4f(self, input_data, is_inverse):
|
||
ifft_ret = self.ifft(input_data)
|
||
if is_inverse:
|
||
ifft_ret = self.fftshift(ifft_ret)
|
||
return ifft_ret
|
||
|
||
PI = 3.1415926535897
|
||
LC = 299.792458
|
||
maichongshu = 400
|
||
jln = 512
|
||
f0 = 9500
|
||
fs = 3600.0
|
||
B = 3000
|
||
PRF = 2000
|
||
vdengxiao = 0
|
||
mubiaojuli = 700e3
|
||
mohuhalfshu = 1
|
||
mohuhalfshu = (int)((15 * 2 / (LC / f0) / PRF) * 1.5) + 1
|
||
mohuhalfshu = 5
|
||
|
||
rfftn = Get2intm(jln)
|
||
afftn = Get2intm(maichongshu)
|
||
datar, datai = read_hrrp_data("hrrp_400_512.dat", maichongshu, jln, rfftn)
|
||
input_data = np.zeros(datar.shape, dtype=np.complex64)
|
||
input_data = datar + 1j * datai
|
||
input_tensor = ms.Tensor(input_data)
|
||
print(input_data.shape)
|
||
print(input_tensor)
|
||
model = Hrrp(maichongshu, jln, f0, fs, B, mohuhalfshu, PRF, vdengxiao, mubiaojuli)
|
||
out = model(input_tensor)
|
||
out = out.asnumpy()
|
||
dr_out = out.real
|
||
di_out = out.imag
|
||
write_hrrp_binary_v2("hrrp_out_400_512_py.dat", dr_out, di_out, maichongshu, jln)
|
||
print(out.shape)
|
||
print(out)
|
||
ms.export(model, input_tensor, file_name="hrrp_400_512", file_format='MINDIR')
|
||
``` |