CN101833468B - 在高性能计算系统中生成向量处理指令集结构的方法 - Google Patents

Abstract

本发明公开了一种在高性能计算系统中生成向量处理指令集结构的方法，该方法包括：步骤1：分析数字信号处理领域内的核心算法，得到各核心算法的逻辑表达式；步骤2：从该逻辑表达式中抽象出数字信号处理领域内的向量元操作；步骤3：根据抽象出的向量元操作，生成针对专用向量指令语法和行为的向量处理指令集结构。利用本发明，生成了适合于特定领域的向量处理指令集，大大提高了处理器的性能。

Description

在高性能计算系统中生成向量处理指令集结构的方法

技术领域

本发明涉及数字信号处理技术领域，特别是一种在高性能计算系统中生成向量处理指令集结构的方法。

背景技术

在过去的四十年期间，DSP处理器得到了突飞猛进的发展。但DSP计算能力必须进一步提升才能满足当前应用需求。半导体工艺是DSP发展的基础，一直是DSP发展的最主要动力，过去DSP性能提升绝大部分来源于芯片工作频率的提高。但由于集成电路制造工艺和功耗的限制，主频提升的时代已经结束，必须依靠架构创新才能进一步提高处理器性能，但架构的优化只能针对特定应用领域。

数字信号处理中绝大部分算法都可以抽象成向量之间的运算。但当今DSP处理器大都以标量为基本处理单位，只对标量计算进行优化设计，向量数据也转化成标量数据后再进行计算。由于工艺和功耗的限制，处理器时钟频率无法进一步提高；当今的处理器架构再也无法通过提高标量数据计算速度来提高向量计算速度。处理器必须针对向量数据的内在并行性进行优化设计。

“指令”是设计者所定义的、芯片可以理解的符号。通过向芯片发送不同的指令序列，程序员指定芯片不同时刻的动作。芯片所能理解的所有指令的集合，即为该芯片的指令集。程序员利用指令集中的指令，实现信号处理算法。

已有一些专利和文章讨论了向量处理器的指令集架构，但这些指令集都只讨论了向量的四则运算，没有针对数字信号处理算法进行分析和提取，所得到的指令集架构不一定适合指定的应用领域；同时，这些专利都是在标量的基础上来讨论向量计算，利用标量的计算资源来实现向量计算，无法进一步利用向量数据内在并行性进行优化设计。

提取和设计向量指令集是进行向量化DSP处理器设计的第一个步骤，指令集结构是否适合应用领域特性将直接影响芯片的整体性能，因此，向量指令集结构设计在芯片设计流程中是一个非常重要的环节。如何将数字信号处理算法抽象成以向量为基本单位的、并行的元操作，同时定义相应的向量指令集，是数字信号处理器研究领域中的一个重要问题。

发明内容

(一)要解决的技术问题

有鉴于此，本发明的主要目的是提供一种在高性能计算系统中生成向量处理指令集结构的方法，以生成适合于特定领域的向量处理指令集，进而提高处理器性能。

(二)技术方案

为达到上述目的，本发明采用的技术方案如下：

一种在高性能计算系统中生成向量处理指令集结构的方法，该方法包括：

步骤1：分析数字信号处理领域内的核心算法，得到各核心算法的逻辑表达式；

步骤2：从该逻辑表达式中抽象出数字信号处理领域内的向量元操作；

步骤3：根据抽象出的向量元操作，生成针对专用向量指令语法和行为的向量处理指令集结构。

上述方案中，步骤1中所述核心算法至少包括脉冲压缩、多普勒处理、空时数据向量生成和三角线性方程组求解。

上述方案中，其特征在于：

对于脉冲压缩，其逻辑表达式为：

$x_1(l,p,n)=\underset{k_p=0}{\overset{K_p-1}{\mathrm{\Σ}}}{h}_p\left(k_p\right)x(l,p,n-k_p)$

其中x(l，p，n)为输入的三维数据，x₁(l，p，n)为输出的三维数据，h_p(k_p)为数值固定的滤波器系数，K_p也为一固定值，表示滤波器系数长度；

对于多普勒处理，其逻辑表达式为：

$x_2(l,k,r)=\underset{p=0}{\overset{P-1}{\mathrm{\Σ}}}d\left(p\right)x_1(l,p,r)e^{j(2\mathrm{\π}/K)\mathrm{pk}}$

其中x₁(l，p，r)为三维输入数据，x₂(l，k，r)为三维输出数据，d(p)为数值固定的窗函数，P也为一固定值，表示窗函数长度；

对于空时数据向量生成算法，其逻辑表达式为：

$\stackrel{\→}{x}(k,r)={[x_2(0,k_{\min },r)...x_2(L-1,k_{\min },r)...x_2(0,k_{\max },r)...x_2(L-1,k_{\max },r)]}^T$

其中x₂(l，k，r)为L×K×R的三维数据，而

为L×Q的二维向量，Q为预定值，k_min＝mod_k(k-(Q-1)/2)，k_max＝mod_k(k+(Q-1)/2)之间，mod_k()表示按k求模；

对于QR分解算法，其逻辑表达式为：

$\stackrel{\→}{A}=\stackrel{\→}{Q}\stackrel{\→}{R}$

其中

为已知矩阵，维度为N×L，且N≥L；

和

为待求解的矩阵，为N×N的单位正交矩阵，而

为N×L维的上三角矩阵；

对于三角线性方程组求解，其逻辑表达式为：

$\stackrel{\→}{R}\stackrel{\→}{X}=\stackrel{\→}{B}$

其中为已知的L×L维上三角矩阵或下三角矩阵，为已知1×L维向量，

为待求解的1×L向量。

上述方案中，所述步骤2包括：根据得到的逻辑表达式，将该核心算法分解为标量操作和向量元操作的组合，进而抽象出数字信号处理领域内的向量元操作。

上述方案中，在将该核心算法分解为标量操作和向量元操作的组合时，首先将该核心算法待处理的对象分解成各个向量单元，然后再尽量将核心算法映射成向量单元之间的元操作。

上述方案中，其特征在于：

对于逻辑表达式为的核心算法，抽象出向量元操作的过程首先对输入数据分段，每段大小与滤波器系数相当，然后用FFT实现每段输入数据的卷积运算，最后通过选择性相加，得到最终结果；在重叠保留法中，将输入数据分成多个互相重叠的段，在最后结果中舍弃重叠部分，组成输出结果；由此可以得出该过程中的向量元操作为：向量FFT变换、向量对应元素相乘、向量IFFT变换；

对于逻辑表达式为

的核心算法，抽象出向量元操作的过程为以d(p)向量为基本长度单位，首先进对输入向量对应元素进行相乘，再对乘法结果进行FFT变换；由此可以得出该过程中的向量元素操作为：向量对应元素相乘，向量FFT变换；

对于逻辑表达式为

的核心算法，抽象出向量元操作的过程为在三维数据输入数据x₂(l，k，r)中，先读取长度为L的向量[x₂(l，k_min，r)]，(l＝0…L-1，r＝0…R-1)，写入目的地址，再读取长度为L的向量[x₂(l，k_min+1，r)]，(l＝0…L-1，r＝0…R-1)，写入目的地址，直到读取到向量[x₂(l，k_max+1，r)]，(l＝0…L-1，r＝0…R-1)，并将其写入目的地址为止；因此，该过程向量元操作为向量离散加载，向量连续存储；

对于逻辑表达式为

的核心算法，抽象出向量元操作是Gram-Schmit正交化，假设m×n的矩阵

包含n个列向量，

设

表示Q矩阵的第i列，R_ij表示R矩阵中的第(i，j)元素，可按如下步骤获得

的分解结果；

首先对向量

进行标准正交化，同时，从中减去与

平行的分量，即：

$\left.\begin{array}{c}{R}_{11}=\left|\right|\stackrel{\→}{a_1}\left|\right|,\stackrel{\→}{q_1}=\stackrel{\→}{a_1}/R_{11}\\ R_{1j}={\stackrel{\→}{q_1}}^H\stackrel{\→}{a_j},{\stackrel{\→}{a_j}}^{\left(1\right)}=\stackrel{\→}{a_j}-\stackrel{\→}{q_1}/R_{1j},2\≤j\≤n\end{array}\right\}$

经过上述运算后，与

正交；然后将

正交化，并从减去与

并行的分量，如下所示：

$\left.\begin{array}{c}{R}_{22}=\left|\right|{\stackrel{\→}{a_2}}^{\left(1\right)}\left|\right|,\stackrel{\→}{q_2}=\stackrel{\→}{a_2}/R_{22}\\ R_{2j}={\stackrel{\→}{q_2}}^H{\stackrel{\→}{a_j}}^{\left(1\right)},{\stackrel{\→}{a_j}}^{\left(2\right)}={\stackrel{\→}{a_j}}^{\left(1\right)}-\stackrel{\→}{q_2}/R_{2j},3\≤j\≤n\end{array}\right\}$

这样构造的向量

均与

正交；重复这一过程，即可得到最后的

矩阵和

矩阵；由以上分析可以看出，QR分解中也以向量为基本运算单位，在进行计算时，也需要从存储器中加载/存储向量单元；因此，QR分解可归结为以下几个方面的操作：求解向量范数，即向量与自身点乘后，开根号；向量与标量进行运算，如

向量点乘，如

向量相减，如

对于逻辑表达式为

的核算法，

为下三角矩阵，抽象出向量元操作的过程为先将方程组写成如下形式：

其中R_ij为下三角矩阵元素，x_i为未知数，b_i为常系数，可利用前向替代法来求解该方程组；首先根据第一行求解x₁，然后将x₁代入第二行，得出x₂，…，x_m依次可解出全部未知数；即

$x_m=\frac{(b_m-\underset{i=1}{\overset{m-1}{\mathrm{\Σ}}}{R}_{\mathrm{mixi}})}{R_{\mathrm{mm}}}$

因此，当m较大时，主要操作为向量点乘

当R为上三角矩阵时，抽取方法与此相同，因此，该类核心算法的向量元操作为向量点积。

上述方案中，所述步骤3包括：统计抽象出的向量元操作，生成针对专用向量指令语法和行为的指令，并将这些指令组合构成向量处理指令集结构。

上述方案中，所述向量处理指令集结构采用VLIW指令结构，至少包括如下四类指令：

标量计算类指令，该类指令字用于实现标量的运算，该标量的运算至少包括标量的加、减、乘、除和逻辑运算；

程序控制类指令，该类指令字用于实现程序控制，改变程序执行流程，至少包括跳转指令、子程序控制指令和循环指令；

向量计算类指令，该类指令用于实现向量运算，至少包括向量对应元素加、减，向量对应元素相乘，向量点积，向量横向求和；

数据搬运类指令，该类指令用于实现标量数据和向量数据的搬运，包括寄存器与寄存器之间、寄存器与存储器之间的数据搬运。

上述方案中，所述标量计算类指令、程序控制类指令和向量计算类指令是单独执行或者并行执行的。

上述方案中，所述标量计算类指令、程序控制类指令、向量计算类指令和数据搬运类指令均由多个指令行构成，每个指令行包括1至4个指令字，一个指令行中的所有指令字并行执行，每个指令字是最终指令集中的任何一条指令，但每个指令行中的所有指令字最多只能对应于一个标量计算类指令、程序控制类指令、向量计算类指令或数据搬运类指令。

(三)有益效果

从上述技术方案可以看出，本发明具有以下有益效果：

1、本发明提供的这种在高性能计算系统中生成向量处理指令集结构的方法，生成的向量处理指令集结构以向量为基本运算对象，能更简洁的表达程序员的意图；向量指令显示表达数据之间的并行，硬件可以更加高效的实现数据并行计算，而不用考虑数据之间的相关；由于向量处理指令以向量为粒度并且向量内元素互不相关，硬件可以进行更加高效的设计，如深度流水、更多并行的功能单元以及更低的时钟频率，因此，生成了适合于特定领域的向量处理指令集，大大提高了处理器的性能。

2、本发明提供的这种在高性能计算系统中生成向量处理指令集结构的方法，从指定领域的算法分析入手，抽取出向量元操作，所得到家的指令非常适合应用领域特性。利用该指令集，可提高应用领域程序执行效率。

3、本发明提供的这种在高性能计算系统中生成向量处理指令集结构的方法，提供了一种基本的VLIW指令架构，在得出“向量相关指令”集合后，可快速生成完备的指令集结构。

附图说明

图1是本发明提供的在高性能计算系统中生成向量处理指令集结构的方法流程图；

图2为本发明提供的在高性能计算系统中生成向量处理指令集结构的示意图；

图3是本发明实施例中相控阵雷达CPI数据的示意图；

图4是本发明实施例中元素空间后多普勒处理自适应算法数据流示意图；

图5是本发明实施例中的组成示意图；

图6是本发明实施例中重叠保留法实现快速卷积的示意图；

图7是本实施例中QR分解过程的示意图；

图8是本实施例中权重相乘的示意图；

图9是本实施例中目标检测算法中采样点定义的示意图；

图10是本实施例中求解极值点的示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

如图1所示，图1是本发明提供的在高性能计算系统中生成向量处理指令集结构的方法流程图，该方法包括以下步骤：

步骤1：分析数字信号处理领域内的核心算法，得到各核心算法的逻辑表达式。

在本步骤中，首先分析当前数字信号处理领域内的核心算法，概括各个核心算法的逻辑表达式，从算法层次了解应用领域的计算需求。核心算法至少包括脉冲压缩、多普勒处理、空时数据向量生成和三角线性方程组求解。

对于脉冲压缩，其逻辑表达式为

其中x(l，p，n)为输入的三维数据，x₁(l，p，n)为输出的三维数据，h_p(k_p)为数值固定的滤波器系数，K_p也为一固定值，表示滤波器系数长度；

对于多普勒处理，其逻辑表达式为：

其中x₁(l，p，r)为三维输入数据，x₂(l，k，r)为三维输出数据，d(p)为数值固定的窗函数，P也为一固定值，表示窗函数长度；

对于空时数据向量生成算法，其逻辑表达式为：

其中x₂(l，k，r)为L×K×R的三维数据，而

为L×Q的二维向量，Q为预定值，k_min＝mod_k(k-(Q-1)/2)，k_max＝mod_k(k+(Q-1)/2)之间，mod_k()表示按k求模；

对于QR分解算法，其逻辑表达式为

其中为已知矩阵，维度为N×L，且N≥L；

和为待求解的矩阵，为N×N的单位正交矩阵，而

为N×L维的上三角矩阵；

对于三角线性方程组求解，其逻辑表达式为其中为已知的L×L维上三角矩阵或下三角矩阵，

为已知1×L维向量，为待求解的1×L向量。

步骤2：从该逻辑表达式中抽象出数字信号处理领域内的向量元操作。

在本步骤中，是根据得到的逻辑表达式，考虑如何在计算机中以向量的方式实现核心算法，将该核心算法分解为标量操作和向量元操作的组合，进而抽象出数字信号处理领域内的向量元操作。在将该核心算法分解为标量操作和向量元操作的组合时，首先将该核心算法待处理的对象分解成各个向量单元，然后再尽量将核心算法映射成向量单元之间的元操作。

对于逻辑表达式为的核心算法，抽象出向量元操作的过程首先对输入数据分段，每段大小与滤波器系数相当，然后用FFT实现每段输入数据的卷积运算，最后通过选择性相加，得到最终结果；在重叠保留法中，将输入数据分成多个互相重叠的段，在最后结果中舍弃重叠部分，组成输出结果；由此可以得出该过程中的向量元操作为：向量FFT变换、向量对应元素相乘、向量IFFT变换；

对于逻辑表达式为

的核心算法，抽象出向量元操作的过程为以d(p)向量为基本长度单位，首先进对输入向量对应元素进行相乘，再对乘法结果进行FFT变换；由此可以得出该过程中的向量元素操作为：向量对应元素相乘，向量FFT变换；

对于逻辑表达式为

的核心算法，抽象出向量元操作的过程为在三维数据输入数据x₂(l，k，r)中，先读取长度为L的向量[x₂(l，k_min，r)]，(l＝0…L-1，r＝0…R-1)，写入目的地址，再读取长度为L的向量[x₂(l，k_min+1，r)]，(l＝0…L-1，r＝0…R-1)，写入目的地址，直到读取到向量[x₂(l，k_max+1，r)]，(l＝0…L-1，r＝0…R-1)，并将其写入目的地址为止；因此，该过程向量元操作为向量离散加载，向量连续存储；

对于逻辑表达式为

的核心算法，抽象出向量元操作是Gram-Schmit正交化，假设m×n的矩阵包含n个列向量，设

表示Q矩阵的第i列，R_ij表示R矩阵中的第(i，j)元素，可按如下步骤获得

的分解结果；

首先对向量

进行标准正交化，同时，从

中减去与

平行的分量，即：

$\left.\begin{array}{c}{R}_{11}=\left|\right|\stackrel{\→}{a_1}\left|\right|,\stackrel{\→}{q_1}=\stackrel{\→}{a_1}/R_{11}\\ R_{1j}={\stackrel{\→}{q_1}}^H\stackrel{\→}{a_j},{\stackrel{\→}{a_j}}^{\left(1\right)}=\stackrel{\→}{a_j}-\stackrel{\→}{q_1}/R_{1j},2\≤j\≤n\end{array}\right\}$

经过上述运算后，

与

正交；然后将

正交化，并从

减去与

并行的分量，如下所示：

$\left.\begin{array}{c}{R}_{22}=\left|\right|{\stackrel{\→}{a_2}}^{\left(1\right)}\left|\right|,\stackrel{\→}{q_2}=\stackrel{\→}{a_2}/R_{22}\\ R_{2j}={\stackrel{\→}{q_2}}^H{\stackrel{\→}{a_j}}^{\left(1\right)},{\stackrel{\→}{a_j}}^{\left(2\right)}={\stackrel{\→}{a_j}}^{\left(1\right)}-\stackrel{\→}{q_2}/R_{2j},3\≤j\≤n\end{array}\right\}$

这样构造的向量

均与

正交；重复这一过程，即可得到最后的

矩阵和

矩阵；由以上分析可以看出，QR分解中也以向量为基本运算单位，在进行计算时，也需要从存储器中加载/存储向量单元；因此，QR分解可归结为以下几个方面的操作：求解向量范数，即向量与自身点乘后，开根号；向量与标量进行运算，如

向量点乘，如

向量相减，如

对于逻辑表达式为

的核心算法，

为下三角矩阵，抽象出向量元操作的过程为先将方程组写成如下形式：

其中R_ij为下三角矩阵元素，x_i为未知数，b_i为常系数，可利用前向替代法来求解该方程组；首先根据第一行求解x₁，然后将x₁代入第二行，得出x₂，…，x_m依次可解出全部未知数；即

$x_m=\frac{(b_m-\underset{i=1}{\overset{m-1}{\mathrm{\Σ}}}{R}_{\mathrm{mixi}})}{R_{\mathrm{mm}}}$

因此，当m较大时，主要操作为向量点乘

当R为上三角矩阵时，抽取方法与此相同，因此，该类核心算法的向量元操作为向量点积。

当然一个完整的算法不可能完全映射成纯粹的向量之间的操作，但如果数字信号处理领域中的某一特定领域内只有很小一部分运算可以映射成向量元操作，而需要很大一部分标量运算，那么该特定领域是不适合进行向量计算的，也无法生成适合该特定领域的向量指令集结构。

步骤3：根据抽象出的向量元操作，生成针对专用向量指令语法和行为的向量处理指令集结构。

在本步骤中，统计抽象出的向量元操作，生成针对专用向量指令语法和行为的指令，并将这些指令组合构成向量处理指令集结构。

但指令集对于数字信号处理领域中的某一特定领域必须是完备的，即该特定领域内的所有算法必须都能映射成指令集中指令的组合。“向量相关指令”的集合只定义了向量处理指令，不满足指令完备性要求。因此，最终的指令集中还必须包括标量计算指令、标量数据搬运指令和程序控制指令。

标量计算、程序控制以及向量计算常常可以并行执行，因此，本发明提供一个最终的指令集结构框架。所述向量处理指令集结构采用VLIW指令结构，至少包括如下四类指令：

1)、标量计算类指令，该类指令字用于实现标量的运算，该标量的运算至少包括标量的加、减、乘、除和逻辑运算等；

2)、程序控制类指令，该类指令字用于实现程序控制，改变程序执行流程，至少包括跳转指令、子程序控制指令和循环指令等；

3)、向量计算类指令，该类指令用于实现向量运算，至少包括向量对应元素加、减，向量对应元素相乘，向量点积，向量横向求和；

4)、数据搬运类指令，该类指令用于实现标量数据和向量数据的搬运，包括寄存器与寄存器之间、寄存器与存储器之间的数据搬运。

最终指令集结构如图2所示。向量相关指令集合进一步归类成“向量计算类指令”和“数据搬运类指令”。标量计算指令、程序控制类指令、向量计算类指令、数据搬运类指令构成“指令集”。标量计算类指令、程序控制类指令、向量计算类指令和数据搬运类指令均由多个指令行构成，每个指令行包括1至4个指令字，一个指令行中的所有指令字并行执行，每个指令字是最终指令集中的任何一条指令，但每个指令行中的所有指令字最多只能对应于一个标量计算类指令、程序控制类指令、向量计算类指令或数据搬运类指令。

下面以相控制雷达信号处理为例，详细说明三个步骤的工作内容。

步骤一：核心算法分析。

相控阵雷达以CPI(Coherent Processing Interval)为数据处理单元，一个CPI数据为三维矩阵，如图3所示，图3是本发明实施例中相控阵雷达CPI数据的示意图。

以x(l，p，n)表示CPI中的一个数据，其中1＝0，…，L-1为通道序号，p＝0，…，P-1，为脉冲序号，n＝0，…，N-1，为回波序号或距离序号。对于大多数机载雷达系统来说，由于计算量和训练数据集太大，无法实现STAP中的完全自适应的最优权重计算。因此，实际代码中都采用部分自适应算法，得到系统的次优权重。以元素空间后多普勒自适应算法为例，算法的主要流程图4所示，图4是本发明实施例中元素空间后多普勒处理自适应算法数据流示意图。Channel表示通道维度，Pulse表示脉冲序号维度，Time表示回波序号维度，Doppler表示速度维度。

下面分别介绍各个算法细节。

1)脉冲压缩

脉冲压缩实际上是一个匹配滤波的过程，用一般的FIR即可实现，假设滤波器系数为h_p(k)，k＝0，…，K_p-1，则滤波器的输出为

$x_1(l,p,n)=\underset{k_p=0}{\overset{K_p-1}{\mathrm{\Σ}}}{h}_p\left(k_p\right)x(l,p,n-k_p)---\left(1\right)$

其中l＝0，…，L-1为通道序号，p＝0，…，P-1，为脉冲序号，n＝0，…，N-1，为回波序号。在实际的代码中，通常采用快速卷积算法实现公式(1)。

2)多普勒处理

在脉冲压缩之后，即多普勒处理。通过多普勒处理，可以得到某一距离上物体的多普勒频率，即物体的相对运动速度。对于某一通道中某一距离点数据，对P个脉冲回波数据进行K阶傅立叶变换得到多普勒频率。同时，在进行变换之前需要对数据进行加窗处理。多普勒处理如下所示：

$x_2(l,p,k)=\underset{p=0}{\overset{P-1}{\mathrm{\Σ}}}d\left(p\right)x_1(l,p,r)e^{j(2\mathrm{\π}/K)\mathrm{pk}}---\left(2\right)$

其中1＝0，…，L-1为通道序号，k＝0，…，K-1，为多普勒频率序号，r＝0，…，N-1，为距离序号，d(p)为加窗函数

3)自适应处理

基于多普勒频率的高阶自适应算法同时利用时间和空间上的自由度进行自适应滤波，得出最优权重。空间上的自由度即L个阵列通道，时间上的自由度即指定的Q个相邻多普勒频率点。

某个距离r和某个多普勒频率点k的权重由时空数据向量

的二次统计结果计算得来。时空数据向量

由L个通道的Q个相邻多普勒频率点数据组成。Q个相邻多普勒频率点的索引值在k_min＝mod_k(k-(Q-1)/2)到k_max＝mod_k(k+(Q-1)/2)之间。mod_k()表示按k求模。因此，

向量长度为L’＝LQ。如下公式表示了的数据组成。

$\stackrel{\→}{x}(k,r)={[x_2(0,k_{\min },r)...x_2(L-1,k_{\min },r)...x_2(0,k_{\max },r)...x_2(L-1,k_{\max },r)]}^T---\left(3\right)$

其中x₂(l，k，r)为多普勒处理后的第1个通道数据的第r个距离单元的第k个多普勒频率点。图5表示了

的组成示意图：

在定义之后，可定义相应的自适应权重向量

$\stackrel{\→}{w}(k,r)={[w(0,k_{\min },r)...w(L-1,k_{\min },r)...w(0,k_{\max },r)...w(L-1,k_{\max },r)]}^T---\left(4\right)$

其中w(l，k，r)表示第1个通道数据的第r个距离单元的第k个多普勒频率点的权重。权重的具体数值需要求解方程，这将在后面讨论。经过自适应处理后的输出即为

$x_3(k,r)={\stackrel{\→}{w}}^H(k,r)\stackrel{\→}{x}(k,r)---\left(5\right)$

其中k＝0，…，K-1，为多普勒频率序号，r＝0，…，N-1，为距离序号。x₃(k，r)构成了最终的速度-距离图，各种目标检测算法即可应用于该数据。

求解权重值之前定义时空向量数据的协方差矩阵为

$\stackrel{\→}{\mathrm{\ψ}}(k,r)=E\left\{\stackrel{\→}{x}(k,r){\stackrel{\→}{x}}^H(k,r)\right\}---\left(6\right)$

其中E{·}表示求数学期望值，则自适应权重可以通过解如下线性方程组获得：

$\stackrel{\→}{\mathrm{\ψ}}(k,r)\stackrel{\→}{w}(k,r)=\mathrm{\γ}\stackrel{\→}{S}---\left(7\right)$

其中γ为缩放因子，使得

为单位阵。而为L’×1的目标时空方向向量。

$\stackrel{\→}{S}={[s\left(\mathrm{0,0}\right)...s(L-\mathrm{1,0})...s(0,Q-1)...s(L-1,Q-1)]}^T---\left(8\right)$

其中

s(l，q)＝s’(l)g(q+mod_k(-(Q-1)/2))    (9)

q＝0，…，Q-1，l＝0，…，L-1，s’(l)为第l个通道中与目标对应的空间方向因子，而

$g\left(k\right)=\underset{p=0}{\overset{P-1}{\mathrm{\Σ}}}d\left(p\right)e^{j(2\mathrm{\π}/K)\mathrm{pk}}---\left(10\right)$

其中d(·)为多普勒处理中FFT窗函数。

从公式    (7)可以看出，该STAP算法依赖于时空协方差矩阵，而实际应用中，通常无法预知该矩阵的值，只能通过估计获得。一般取相邻距离单元的协方差平均值。将N个距离单元分成M个组，每组中包含N_R个连续的距离单元，即M＝N/N_R。对时空向量

的内积取平均，即为第k个多普勒频率点的第m个距离组的协方差矩阵.

$\underset{\‾}{\overset{\→}{\mathrm{\ψ}}}(k,m)=\frac{1}{N_R}\underset{r=r_1}{\overset{r_1+N_R-1}{\mathrm{\Σ}}}\stackrel{\→}{x}(k,r){\stackrel{\→}{x}}^H(k,r)---\left(11\right)$

其中

为协方差矩阵的估计值，r₁＝mN_R，m＝0，…，M-1，而k＝0，…，K-1。协方差矩阵的估计值将代替公式(7)中的协方差矩阵，用于计算自适应权重向量。第m个组中所有距离单元的第k个多普勒频率点都对应该权重向量。

第m个组中的N_R个时空向量构成一个时空矩阵可以得到更简洁的表达为：

$\stackrel{\→}{X}(k,m)=[\stackrel{\→}{x}(k,m{N}_R)...\stackrel{\→}{x}(k,(m+1)N_R-1)]---\left(12\right)$

$\underset{\‾}{\overset{\→}{\mathrm{\ψ}}}(k,m)=\frac{1}{N_R}\stackrel{\→}{X}(k,m){\stackrel{\→}{X}}^H(k,m)---\left(13\right)$

则公式    (7)变为

$\stackrel{\→}{X}(k,m){\stackrel{\→}{X}}^H(k,m)\stackrel{\→}{w}(k,m)=\mathrm{\γ}{N}_R\stackrel{\→}{S}---\left(14\right)$

其中

为第m个组中所有距离单元的第k个多普勒频率点的权重向量。为了求解公式(14)公式中的

首先对L’×N_R维的矩阵进行QR分解，使得

其中

为N_R×N_R的单位正交矩阵，而为N_R×L’维的上三角矩阵。通常要求N_R＞L’，从而

可以写成

的形式，

为L’×L’维的满秩上三角矩阵。对

进行QR分解后，则有

$\stackrel{\→}{X}(k,m){\stackrel{\→}{X}}^H(k,m)={\stackrel{\→}{R}}^T{\stackrel{\→}{Q}}^T{\stackrel{\→}{Q}}^{*}{\stackrel{\→}{R}}^{*}={{\stackrel{\→}{R}}_1}^T{{\stackrel{\→}{R}}_1}^{*}---\left(15\right)$

将公式(15)代入(14)可得到

${{\stackrel{\→}{R}}_1}^T{{\stackrel{\→}{R}}_1}^{*}\stackrel{\→}{w}(k,m)=\mathrm{\γ}{N}_R\stackrel{\→}{S}---\left(16\right)$

令

则公式(16)变为

${{\stackrel{\→}{R}}_1}^T\stackrel{\→}{p}=\mathrm{\γ}{N}_R\stackrel{\→}{S}---\left(17\right)$

利用前向消除法求解该线性方程组可得向量，再利用后向替代法求解线程方程组

${{\stackrel{\→}{R}}_1}^{*}\stackrel{\→}{w}(k,m)=\stackrel{\→}{p}---\left(18\right)$

可得到最终的权重向量

在得到权重向量后，将其与第m个距离组的第k个多普勒频率点数据相乘，得到N_R个最终输出结果。如下以公式所示：

$x_4(k,m)={\stackrel{\→}{w}}^H(k,m)\stackrel{\→}{X}(k,m)---\left(19\right)$

对M个距离组的K个多普勒频率点数据重复这一过程，可得到最终的多普勒-距离图。所以一共需要M×K个求解过程才能得到最终结果。

步骤二：向量元操作抽取。雷达信号处理对象为CPI立方体数据，操作步骤都针对立方体数据中一行或一列，因此，处理对象可方便的分解成行向量或列向量。在进行计算时，都需要先从存储器中读取向量单元，计算完成后，再将计算结果写回存储器。同时，由于CPI为多维数据，在存储器中只能按行存储或按列存储。当按行存储按行读取时，向量的各个元素地址是连续的。但如果按列存储按行读取时，各个向量之间的地址是不连续的。因此，各个计算步骤都需要进行如下四个数据搬运操作：

向量连续地址加载

向量离散地址加载

向量连续地址存储

向量离散地址存储

1)脉冲压缩

脉冲压缩即匹配滤波过程，在算法上等同于FIR滤波，通过FFT变换后在频率上进行滤波可减少算法运算量。但脉冲压缩中滤波器系数长度要比输入数据长度小得多，在FFT运算之前需要对滤波器系数补0，反而影响算法性能。因此，实际应用中都采用快速卷积实现FIR脉冲压缩。快速卷积即首先对输入数据分段，每段大小与滤波器系数相当，然后用FFT实现每段输入数据的卷积运算，最后通过选择性相加，得到最终结果。在重叠保留法中，将输入数据分成多个互相重叠的段，如图6中的x₁(n)，x₂(n)等所示，在最后结果中舍弃重叠部分，组成输出结果。图6中H(k)为匹配滤波器系数的FFT变换，可事先计算好存放在内存中，X(k)为各个向量段x_j(n)的FFT变换，Y(k)为X(k)和H(k)对应元素相乘后的结果，y_j(n)为各个Y(k)的FFT逆变换。

可以看出，快速卷积中以向量为基本单位，进行以下几个方面的操作：

向量FFT变换

向量对应元素相乘

向量IFFT变换

2)多普勒处理

而在多普勒处理中，也以向量为基本单位，基本操作为加窗操作和FFT变换，而加窗操作即向量对应元素相乘。因此，该步骤可以归结为：

向量FFT变换

向量对应元素相乘

3)自适应处理

自适应处理主要包括QR分解、前向消除求解线性方程组、后向代替求解线性方程组、权重相乘四个小步骤。

QR分解主要操作是Gram-Schmit正交化，假设m×n的矩阵

包含n个列向量，

设

表示Q矩阵的第i列，R_ij表示R矩阵中的第(i，j)元素，可按如下步骤获得

的分解结果。

首先对向量进行标准正交化，同时，从

中减去与平行的分量，即：

$\left.\begin{array}{c}{R}_{11}=\left|\right|\stackrel{\→}{a_1}\left|\right|,\stackrel{\→}{q_1}=\stackrel{\→}{a_1}/R_{11}\\ R_{1j}={\stackrel{\→}{q_1}}^H\stackrel{\→}{a_j},{\stackrel{\→}{a_j}}^{\left(1\right)}=\stackrel{\→}{a_j}-\stackrel{\→}{q_1}/R_{1j},2\≤j\≤n\end{array}\right\}$

经过上述运算后，

与

正交。然后将

正交化，并从

减去与并行的分量，如下所示：

$\left.\begin{array}{c}{R}_{22}=\left|\right|{\stackrel{\→}{a_2}}^{\left(1\right)}\left|\right|,\stackrel{\→}{q_2}=\stackrel{\→}{a_2}/R_{22}\\ R_{2j}={\stackrel{\→}{q_2}}^H{\stackrel{\→}{a_j}}^{\left(1\right)},{\stackrel{\→}{a_j}}^{\left(2\right)}={\stackrel{\→}{a_j}}^{\left(1\right)}-\stackrel{\→}{q_2}/R_{2j},3\≤j\≤n\end{array}\right\}$

这样构造的向量均与

正交。重复这一过程，即可得到最后的

矩阵和

矩阵。如图7所示，图7是本实施例中QR分解过程的示意图。

由以上分析可以看出，QR分解中也以向量为基本运算单位，在进行计算时，也需要从存储器中加载/存储向量单元。因此，QR分解可归结为以下几个方面的操作：

求解向量范数，即向量与自身点乘后，开根号。

向量与标量进行运算，如

向量点乘，如

向量相减，如

在求得

矩阵后，需要求解线程方程组(17)和(18)，由于

为下三角矩阵，可将方程组(17)写成如下形式：

其中R_ij为下三角矩阵元素，x_i为未知数，b_i为常系数，可利用前向替代法来求解该方程组。首先根据第一行求解x₁，然后将x₁代入第二行，得出x₂，…，x_m依次可解出全部未知数。即

$x_m=\frac{(b_m-\underset{i=1}{\overset{m-1}{\mathrm{\Σ}}}{R}_{\mathrm{mixi}})}{R_{\mathrm{mm}}}---\left(20\right)$

因此，当m较大时，主要操作为向量点乘公式(18)求解方法与此相同。因此该步骤主要操作为向量点乘，在求得权重后，需要将权重应用于输入数据，即进行权重相乘。在STAP中，权重相乘即权重向量与输入矩阵相乘，得到距离组中某个多普勒频率点上信号的幅值。如图8所示，图8是本实施例中权重相乘的示意图。因此，主要运算类型为向量点乘。

4)目标检测

恒虚警目标检测算法在距离-多普勒平面上，在待检测点周围选择连续的采样点，求得采样点平均值后，将待检测点信号幅值与平均值相除。因此，最基本的操作即采样点求和操作。在实际应用中，通常将采样点设定为待检测点周围的矩形区域，如图9所示，图9是本实施例中目标检测算法中采样点定义的示意图。

将矩形区域的每条边看成一向量，则目标检测中主要操作类型为向量横向求和操作。

在恒虚警处理之后，需要从距离-多普勒平面上确定哪些点为目标，以及确定目标参数。即搜索距离-多普勒平面上的极值点，该操作可以分成两个步骤，第一步骤即求各个距离点上的最大多普勒频率，得到距离-最大多普载频率值向量T。第二步即求向量T的极值点。如图10所示，图10是本实施例中求解极值点的示意图。

第一个步骤可以看成是求取向量大值过程，为一向量操作。第二个步骤为一标量处理过程。因此，目标检测中主要向量操作类型为：

向量横向求和

求取向量最大值

步骤三：向量相关指令集定义

总结出相控制阵雷达信号处理中向量的元操作，如下表所示：

雷达信号处理步骤	主要向量操作
		脉冲压缩	向量连续地址加载向量连续地址存储向量FFT向量点乘向量IFFT
多普勒处理	向量连续地址加载向量连续地址存储向量离散地址加载向量离散地址存储向量点乘向量FFT
		自适应处理	向量连续地址加载向量连续地址存储向量离散地址加载向量离散地址存储向量点乘向量相减向量与标量运算向量横向求和

目标检测	向量连续地址加载向量连续地址存储向量离散地址加载向量离散地址存储向量点乘向量最大值

由于向量FFT和向量IFFT包括多个步骤的计算，可将其进一步分解为更细的向量元操作。对于一个长度为N的向量FFT算法，一共需要N/2个旋转因子，进行log2(N)级蝶形运算，每级蝶形运算需要N/2个复数乘法、N个复数加法。对于长度为N的FFT算法，每次蝶形运算中可以将输入数据分解成两部分：一半不参与乘法运算，用向量

表示该部分数据；另一半参与乘法运算，用向量表示该部分数据。另设每级的旋转因子向量为

则每级蝶形可表达为：

因而，向量FFT和向量IFFT也转化成如下向量元素操作的组合：

向量对应元素加

向量对应元素减

向量对应元素乘

在指令集结构中，约定标量寄存器以R表示，向量寄存器以VR表示，向量寄存器的元素个数为N，寄存器下标以s或m或n表示。Mem[i]表示存储器中地址为i的数据内容。则“向量相关的指令”集合如下表所示：

指令类型	指令格式	指令行为(C语言表示)
			向量加	VRs＝VRm+VRn	for(i＝0；i＜N；i++)VRs[i]＝VRm[i]+VRn[i]
向量减	VRs＝VRm-VRn	for(i＝0；i＜N；i++)VRs[i]＝VRm[i]-VRn[i]
			向量乘	VRs＝VRm*VRn	for(i＝0；i＜N；i++)VRs[i]＝VRm[i]*VRn[i]

向量横向求和	Rs＝SUM(VRm)	Rs＝0for(i＝0；i＜N；i++)Rs＝Rs+VRm[i]
			向量点积	Rs＝VRm DOT VRn	Rs＝0for(i＝0；i＜N；i++)Rs＝Rs+VRm[i]*VRn[i]
向量最大值	Rs＝MAX(VRm)	Rs＝VRm[0]for(i＝0；i＜N；i++)Rs＝Rs＞VRm[i]？Rs：VRm[i]
			向量最小值	Rs＝MIN(VRm)	Rs＝VRm[0]for(i＝0；i＜N；i++)Rs＝Rs＜VRm[i]？Rs：VRm[i]
向量连续地址加载	VRs＝[Rm]	for(i＝0；i＜N；i++)VRs[i]＝Mem[Rm+i]
			向量连续地址存储	[Rs]＝VRm	for(i＝0；i＜N；i++)Mem[Rs+i]＝VRs[i]
向量离散地址加载	VRs＝[VRm]	for(i＝0；i＜N；i++)VRs[i]＝Mem[VRm[i]]
			向量离散地址存储	[VRs]＝VRm	for(i＝0；i＜N；i++)Mem[VRs[i]]＝VRm[i]

根据该相量相关指令集合，结合预定义的向量指令集框架构，即可得到适合于相控制阵雷达信号处理的向量指令集结构。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种在高性能计算系统中生成向量处理指令集结构的方法，其特征在于，该方法包括： 

步骤1：分析数字信号处理领域内的核心算法，得到各核心算法的逻辑表达式； 

步骤2：从该逻辑表达式中抽象出数字信号处理领域内的向量元操作； 

步骤3：根据抽象出的向量元操作，生成针对专用向量指令语法和行为的向量处理指令集结构； 

其中，所述步骤2包括：根据得到的逻辑表达式，将该核心算法分解为标量操作和向量元操作的组合，进而抽象出数字信号处理领域内的向量元操作；在将该核心算法分解为标量操作和向量元操作的组合时，首先将该核心算法待处理的对象分解成各个向量单元，然后再利用向量单元的操作特征，评估能否将核心算法映射成向量单元之间的元操作，如该操作能用硬件电路实现，则将操作定义成一条指令。 

2.根据权利要求1所述的在高性能计算系统中生成向量处理指令集结构的方法，其特征在于，步骤1中所述核心算法至少包括脉冲压缩、多普勒处理、空时数据向量生成和三角线性方程组求解。 

3.根据权利要求2所述的在高性能计算系统中生成向量处理指令集结构的方法，其特征在于： 

对于脉冲压缩，其逻辑表达式为

其中x(l，p，n)为输入的三维数据，x₁(l，p，n)为输出的三维数据，h_p(k_p)为数值固定的滤波器系数，K_p也为一固定值，表示滤波器系数长度； 

对于多普勒处理，其逻辑表达式为： 

其中x₁(l，p，r)为三维输入数据， x₂(l，k，r)为三维输出数据，d(p)为数值固定的窗函数，P也为一固定值，表示窗函数长度； 

对于空时数据向量生成算法，其逻辑表达式为： 

其中x₂(l，k，r)为L×K×R的三维数据，而为L×Q的二维向量，Q为预定值，k_min＝mod_k(k-(Q-1)/2)，k_max＝mod_k(k+(Q-1)/2)之间，mod_k()表示按k求模； 

对于QR分解算法，其逻辑表达式为其中

为已知矩阵，维度为N×L，且N≥L；

和

为待求解的矩阵，

为N×N维的单位正交矩阵，而

为N×L维的上三角矩阵； 

对于三角线性方程组求解，其逻辑表达式为

其中

为已知的L×L维上三角矩阵或下三角矩阵，

为已知1×L维向量，

为待求解的1×L维向量。 

4.根据权利要求1所述的在高性能计算系统中生成向量处理指令集结构的方法，其特征在于， 

对于逻辑表达式为

的核心算法，抽象出向量元操作的过程首先对输入数据分段，每段大小与滤波器系数相当，然后用FFT实现每段输入数据的卷积运算，最后通过选择性相加，得到最终结果；在重叠保留法中，将输入数据分成多个互相重叠的段，在最后结果中舍弃重叠部分，组成输出结果；由此可以得出该过程中的向量元操作为：向量FFT变换、向量对应元素相乘、向量IFFT变换； 

对于逻辑表达式为

的核心算法，抽象出向量元操作的过程为：以d(p)向量的尺寸为基本长度单位，截取输入向量x₁(l，p，r)，然后将向量d(p)与截取后输入向量的对应元素进行相乘，再 对乘法结果进行FFT变换；由此可以得出该过程中的向量元素操作为：向量对应元素相乘，向量FFT变换； 

对于逻辑表达式为 

的核心算法，抽象出向量元操作的过程为在三维数据输入数据x₂(l，k，r)中，先读取长度为L的向量[x₂(l，k_min，r)]，(l＝0...L-1，r＝0...R-1)，写入目的地址，再读取长度为L的向量[x₂(l，k_min+1，r)]，(l＝0...L-1，r＝0...R-1)，写入目的地址，直到读取到向量[x₂(l，k_max+1，r)]，(l＝0...L-1，r＝0...R-1)，并将其写入目的地址为止；因此，该过程向量元操作为向量离散加载，向量连续存储； 

对于逻辑表达式为

的核心算法，抽象出向量元操作是Gram-Schmit正交化，假设m×n的矩阵包含n个列向量，

设表示Q矩阵的第i列，R_ij表示R矩阵中的第(i，j)元素，可按如下步骤获得

的分解结果； 

首先对向量

进行标准正交化，同时，从

中减去与

平行的分量，即： 

经过上述运算后，与

正交；然后将

正交化，并从

减去与

并行的分量，如下所示： 

这样构造的向量

均与

正交；重复这一过程，即可得到最后的

矩阵和

矩阵；由以上分析可以看出，QR分解中也以向量为基本运算单位，在进行计算时，也需要从存储器中加载/存储向量单 元；因此，QR分解可归结为以下几个方面的操作：求解向量范数，即向量与自身点乘后，开根号；向量与标量进行运算，如向量点乘，如向量相减，如

对于逻辑表达式为

的核心算法，

为下三角矩阵，抽象出向量元操作的过程为先将方程组写成如下形式： 

其中R_ij为下三角矩阵元素，x_i为未知数，b_i为常系数，可利用前向替代法来求解该方程组；首先根据第一行求解x₁，然后将x₁代入第二行，得出x₂，...，x_m依次可解出全部未知数；即 

因此，当m较大时，主要操作为向量点乘

当R为上三角矩阵时，抽取方法与此相同，因此，该类核心算法的向量元操作为向量点积。 

5.根据权利要求1所述的在高性能计算系统中生成向量处理指令集结构的方法，其特征在于，所述步骤3包括： 

统计抽象出的向量元操作，生成针对专用向量指令语法和行为的指令，并将这些指令组合构成向量处理指令集结构。 

6.根据权利要求1或5所述的在高性能计算系统中生成向量处理指令集结构的方法，其特征在于，所述向量处理指令集结构采用VLIW指令结构，至少包括如下四类指令： 

标量计算类指令，该类指令字用于实现标量的运算，该标量的运算至少包括标量的加、减、乘、除和逻辑运算； 

程序控制类指令，该类指令字用于实现程序控制，改变程序执行流程， 至少包括跳转指令、子程序控制指令和循环指令； 

向量计算类指令，该类指令用于实现向量运算，至少包括向量对应元素加、减，向量对应元素相乘，向量点积，向量横向求和； 

数据搬运类指令，该类指令用于实现标量数据和向量数据的搬运，包括寄存器与寄存器之间、寄存器与存储器之间的数据搬运。 

7.根据权利要求6所述的在高性能计算系统中生成向量处理指令集结构的方法，其特征在于，所述标量计算类指令、程序控制类指令和向量计算类指令是单独执行或者并行执行的。 

8.根据权利要求6所述的在高性能计算系统中生成向量处理指令集结构的方法，其特征在于，所述标量计算类指令、程序控制类指令、向量计算类指令和数据搬运类指令均由多个指令行构成，每个指令行包括1至4个指令字，一个指令行中的所有指令字并行执行，每个指令字是最终指令集中的任何一条指令，但每个指令行中的所有指令字最多只能对应于一个标量计算类指令、程序控制类指令、向量计算类指令或数据搬运类指令。 

Images