# HRRP（高分辨率距离像算法） — MindSpore Signal+ 使用手册 alpha 文档

# HRRP（高分辨率距离像算法）

## 算法概述

高分辨率距离像（High Resolution Range Profile, HRRP）是雷达信号处理中的核心技术，通过对雷达回波信号进行脉冲压缩和多普勒处理，获取目标在距离维上的散射点分布特征。HRRP具有以下特点：

- 距离高分辨：能够区分目标上相邻较近的散射中心

- 特征稳定：为目标识别提供物理特征描述

- 计算高效：相比二维成像计算复杂度较低

- 实时性强：适用于实时目标识别系统

## 代码来源

## MindSpore Signal+ 实现

将原有代码重构为基于计算图的MindSpore Cell，实现雷达信号处理流程模块化。

### 1 数据读取

在Python中使用MindSpore Signal+时，我们可以使用NumPy的fromfile接口读二进制格式文件，因此在Python代码开头需要导入NumPy库。

示例代码：

    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

读到的数据需要转换为复数形式，以便后续处理。

    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__`函数。 示例代码：

    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 相位补偿因子计算

    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变换参数预计算

    # 关键参数计算
    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 泰勒窗函数预计算

    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 循环移位索引预计算

    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处理中，我们只关心每个距离单元的最强散射点，而原始数据的第一行被用来存储这个信息。

核心算法整体流程：

    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 相位补偿

    out = input_data * self.phase_comp

物理意义：

- 补偿目标：校正由于雷达与目标相对运动引起的二次相位误差

- 误差来源：目标斜距变化引起的时间延迟变化

- 数学形式：补偿因子为 exp(j·π·Ka·t²)，其中Ka为调频斜率

#### 3.2 多模糊数处理循环

1. 模糊数范围确定：

    mohushu = 2 * self.mohuhalfshu + 1

设计原理：

- 由于多普勒频率模糊，真实的调频斜率有多个候选值

- 搜索范围：[-mohuhalfshu, +mohuhalfshu]

- 典型设置：mohuhalfshu=5 → 搜索11个候选值

2. 循环处理每个模糊数：

    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变换是一种距离徙动校正算法，用于解决以下问题：

变换公式： 原始坐标：t (方位时间), f (距离频率) 变换后：t' = t × (f₀/(f₀+f))

2. 代码实现

    # 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生成

    # 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 熵值评估

    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 最优结果选择

    # 合并所有结果
    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 端部署模板。具体使用方法可参考 HelloDSP MindSpore Lite端。

以下是完整的MindSpore Signal+实现的代码：

    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')