一种中子脉冲序列的频谱分析方法

摘要

本发明涉及一种中子脉冲序列的频谱分析方法，包括以下步骤：(1)对中子脉冲数据进行采集，得到中子源以及被中子源激发的被探测体产生的中子探测计数的时间分布；(2)对采集后的数据包按一定的长度进行分块，得到各块中的中子脉冲原序列；各数据块之间连续无间隔；(3)对各数据块的数据进行相关计算，得到各数据块或块间的相关序列；(4)对各块或块间的相关序列进行进一步计算得到功率谱；其特征是：所述步骤(3)中对各块数据进行相关计算时使用快速偏移的方法。这种频谱分析方法利用可实时计算的快速偏移相关算法以及功率谱计算步骤的优化，使得该方法能实现中子脉冲序列的实时频谱分析。

说明书

技术领域

本发明涉及一种中子脉冲序列频谱分析方法。

背景技术

中子源驱动噪声分析测量法，也称功率谱密度分析法(power spectral density analysismeasurements)，它的基本原理是：首先创立目标参数的计算公式，其次依托实验测得原函数，并快速计算相关函数和功率谱，最后将这些功率谱代入所创立的目标参数公式，由此便可得到中子辐射环境中相应的目标参数。通过构建实验测量系统进行高达1GHz的采样率进行高速数据采集，可以获得原函数，即中子脉冲序列。

对中子脉冲序列的频谱进行经典分析的方法包括自相关法和周期图法。自相关法是建立在维纳-幸钦定理的基础上的，其基本步骤是先对随机信号的N点数据作自相关估计然后对其作傅立叶变换，便得到功率谱估计。其中，相关函数的计算有两种方法，一是利用公式直接计算，二是通过FFT来实现快速计算。周期图法是先求N点数据的离散时间傅立叶变换，得到X(e^jω)，然后取其幅频特性的平方乘以1/N作为功率谱估计，称为周期图。但是无论是自相关法还是周期图法的计算速度都很慢，不能满足实时分析的要求。例如：当采样率为1GHz，计数率为3×10⁶S^-1(即每秒有3×10⁶个脉冲出现)，块(Block)长度为1024，采集1×10⁹个Block时，要求频谱分析在1小时内完成，也即单个Block的计算分析必须在3.6μs内完成。而在现有PC机平台(Pentium D3.4G，单线程)下，对于单个Block(长度为1024)而言，采用目前使用较为广泛且被公认为最快速的快速傅立叶变换(FFT)程序包之一的FFTW的计算时间为13.5μs，而自相关计算时间则需要130μs。显然，传统的分析方法完全无法满足系统在3.6μs内实时地进行分析的要求。

发明内容

本发明的目的在于提供一种能够对中子脉冲序列进行实时分析的方法。

本发明的目的是通过以下技术方案实现的：一种中子脉冲序列的频谱分析方法，包括以下步骤：(1)对中子脉冲数据进行采集，得到中子源以及被中子源激发的被探测体产生的中子探测计数的时间分布，所述时间分布的形式是由“0”或“1”组成的中子脉冲序列，所述中子脉冲序列按脉冲出现的位置进行存储，没有脉冲的不存储；(2)对采集后的数据包按一定的长度(设为n)进行分块(共m块)，得到各块中的中子脉冲原序列；各数据块之间连续无间隔；(3)对各数据块的数据进行相关计算，得到各数据块或块间的相关序列；(4)对各块或块间的相关序列进行进一步计算得到功率谱；其特征是：所述步骤(3)中对各块数据进行相关计算时使用快速偏移的方法。

若需要得到自相关功率谱，则上述步骤(3)中的相关计算为自相关，所述快速偏移的方法包括以下步骤：a、将上述步骤(2)中所得的中子脉冲原序列x(i)，i＝0，1，...n-1复制一份，得到偏移序列y(i)，y(i)＝x(i)，i＝0，1，L，n-1；b、计算偏移0时中子脉冲原序列x(i)与偏移序列y(i)的相关值，即r(0)＝p，p为中子脉冲原序列x(i)中脉冲为“1”的个数；c、计算中子脉冲原序列x(i)与偏移序列y(i)中数据“1”间的间距，得到最小间距d(d>0)；设偏移计数为k，k＝0；d、计算偏移序列y(i)偏移最小间距d时与中子脉冲原序列x(i)的相关值，即$r\left(d\right)=\underset{i=0}{\overset{n-1-k-d}{\Sigma }}x\left(i\right)y(i+k+d);$偏移计数k＝k+d；e、如果k<n-1，则计算序列y(i)，i＝k，L，n-1与序列x(i)，i＝0，L，n-1-k中数据“1”间的间距，得到此时的最小间距d，如果d>0，转至步骤d；如果k≥n-1或d＝0，则计算完毕，转至下步；f、未计算的其他相关值为0，可预先初始化为0；g、进行计算结果的处理，得到相关序列

对中子脉冲数据进行采集时，中子源以及被中子源激发的被探测体产生的中子由中子探测器探测，利用置于计算机中的高速数据采集卡对源通道信号和被探测的二个通道共三路信号进行高速数据采集，从而得到所述时间分布。当然，被探测的通道也可以为三个或四个。

若需要得到互相关功率频谱，则上述步骤(3)中的相关计算为互相关计算，所述快速偏移的方法包括以下步骤：a、用其中两个通道的序列x₁(i)与x₂(i)计算偏移0时x₁(i)与x₂(i)的相关值，即$r\left(0\right)=\underset{i=0}{\overset{n-1}{\Sigma }}{x}_1\left(i\right)x_2\left(i\right);$b、计算x₁(i)与x₂(i)中数据“1”间的间距，得到最小间距d(d>0)；设偏移计数为k，k＝0；c、计算x₂(i)偏移最小间距d时与x₁(i)的相关值，即$r\left(d\right)=\underset{i=0}{\overset{n-1-k-d}{\Sigma }}{x}_1\left(i\right)x_2(i+k+d);$同时，累计偏移计数k＝k+d；d、如果k<n-1，则计算x₂(i)，i＝k，L，n-1与x₁(i)，i＝0，L，n-1-k中数据“1”间的间距，得到此时的最小间距d；如果d>0，转至步骤d；如果k≥n-1或d＝0，则计算完毕，转至下步；e、未计算的其他相关值为0，可预先初始化为0；f、改变移位方向，重复步骤c、d、e，可算出另一半结果。由此，可以计算得到完整的互相关序列${\hat{R}}_{x1x2}\left(j\right),j=-n+1,...,-\mathrm{1,0,1},...,n-1.$

为进一步缩短计算时间，所述步骤(4)按照如下步骤进行：a、当各数据块中的数累计达到预先设定的阈值时，对各数据块中计算所得的相关序列进行求和得到b、对求和结果进行快速傅立叶变换从而得到功率谱。

所述高速数据采集卡的采样率可以为1GHz，当然也可以为500MHz或者其他频率。

所述数据块的长度可以为1024个字节，当然，块的长度也可以是256、512、2048等。

本发明的有益效果是：

(1)由于在进行相关计算时，利用了中子脉冲数据计数率低的特点，使用快速偏移的方法，省去了没有必要的偏移和相关计算，因此大大减少了相关计算的时间，从而能够对中子脉冲序列进行实时分析；当采样率为1GHz，计数率为3×106S^-1(即每秒有3×106个脉冲出现)，数据块长度为1024，采集1×109个数据块时，利用现有的PC平台(Pentium D3.4G)，频谱分析能在1小时内完成，也即单个Block的计算分析可以在3.6μs内完成。与相应的英特尔数学内核库(Intel Math Kernel Library，MKL)相比较，快速偏移的自相关计算方法的计算效率提高了249倍，具体测试数据见表1。

表1

(2)由于对中子脉冲序列进行数据分块，对各块数据先做相关计算，然后对各块数据相关计算的结果累计求和，再对累计求和结果进行傅立叶变换得到功率谱；而现有的算法是先对各块数据进行相关计算，然后对各相关计算的结果进行傅立叶变化，再对各傅立叶变换的结果求和，因此本发明所采用的方法大大节约了计算的时间，从而能够对中子脉冲序列进行实时分析。

附图说明

图1为本发明实施例1优化功率谱计算示意图；

图2为基于PC平台的频谱分析结构示意图；

图3为中子脉冲序列的频谱测量与分析流程图；

图4为本发明实施例1测试中自相关计算结果图；

图5为本发明实施例1测试中自功率频谱计算结果图；

图6为本发明实施例2优化功率谱计算示意图；

图7为本发明实施例2测试中互相关计算结果图；

图8为本发明实施例2测试中互功率频谱计算结果图。

具体实施方式

下面通过实施例的方式进一步说明本发明，并不因此将本发明限制在所述的实施例范围之中。

实施例1：如图1至图3所示，一种中子脉冲序列的频谱分析方法，包括以下步骤：(1)对中子脉冲数据进行采集，得到中子源以及被中子源激发的被探测体产生的中子探测计数的时间分布，所述时间分布的形式是由“0”或“1”组成的中子脉冲序列，所述中子脉冲序列按脉冲出现的位置进行存储，没有脉冲的不存储；(2)对采集后的数据包按1024个字节的长度(设为n)进行分块(共m块)，得到各块中的中子脉冲原序列x(i)，i＝0，1，...n-1；各数据块之间连续无间隔；(3)对各数据块的数据进行相关计算，得到各数据块或块间的相关序列；(4)对各块或块间的相关序列进行进一步计算得到功率谱。

所述步骤(3)中对各块数据进行的相关计算为自相关计算，使用快速偏移的方法，包括以下步骤：a、将上述步骤(2)中所得的中子脉冲原序列x(i)，i＝0，1，...n-1复制一份，得到偏移序列y(i)，y(i)＝x(i)，i＝0，1，L，n-1；b、计算偏移0时中子脉冲原序列x(i)与偏移序列y(i)的相关值，即r(0)＝p，p为中子脉冲原序列x(i)中脉冲为“1”的个数；c、计算中子脉冲原序列x(i)与偏移序列y(i)中数据“1”间的间距，得到最小间距d(d>0)；设偏移计数为k，k＝0；d、计算偏移序列y(i)偏移最小间距d时与中子脉冲原序列x(i)的相关值，即$r\left(d\right)=\underset{i=0}{\overset{n-1-k-d}{\Sigma }}x\left(i\right)y(i+k+d);$偏移计数k＝k+d；e、如果k<n-1，则计算序列y(i)，i＝k，L，n-1与序列x(i)，i＝0，L，n-1-k中数据“1”间的间距，得到此时的最小间距d，如果d>0，转至步骤d；如果k≥n-1或d＝0，则计算完毕，转至下步；f、未计算的其他相关值为0，可预先初始化为0；g、进行计算结果的处理，得到相关序列

所述步骤(4)按照如下步骤进行：a、当各数据块中的数累计达到预先设定的阈值(如1×107)时，对各数据块中计算所得的自相关序列进行求和得到b、对求和结果进快速傅立叶变换从而得到功率谱。

对中子脉冲数据进行采集时，中子源以及被中子源激发的被探测体产生的中子由中子探测器探测，利用置于计算机中的高速数据采集卡对源通道信号和被探测的二个通道共三路信号进行高速数据采集，从而得到所述时间分布。所述高速数据采集卡的采样率为1GHz，计数率为3×106S^-1。

上述计算完毕后，显示终端的显示如图4和图5所示。

实施例2：如图6所示，本实施例与实施例1所不同的是，则上述步骤(3)中的相关计算为互相关计算，所述快速偏移的方法包括以下步骤：a、用其中两个通道的序列x₁(i)与x₂(i)计算偏移0时x₁(i)与x₂(i)的相关值，即$r\left(0\right)=\underset{i=0}{\overset{n-1}{\Sigma }}{x}_1\left(i\right)x_2\left(i\right);$b、计算x₁(i)与x₂(i)中数据“1”间的间距，得到最小间距d(d>0)；设偏移计数为k，k＝0；c、计算x₂(i)偏移最小间距d时与x₁(i)的相关值，即$r\left(d\right)=\underset{i=0}{\overset{n-1-k-d}{\Sigma }}{x}_1\left(i\right)x_2(i+k+d);$同时，累计偏移计数k＝k+d；d、如果k<n-1，则计算x₂(i)，i＝k，L，n-1与x₁(i)，i＝0，L，n-1-k中数据“1”间的间距，得到此时的最小间距d；如果d>0，转至步骤d；如果k≥n-1或d＝0，则计算完毕，转至下步；e、未计算的其他相关值为0，可预先初始化为0；f、改变移位方向，重复步骤c、d、e，可算出另一半结果。由此，可以计算得到完整的互相关序列${\hat{R}}_{x1x2}\left(j\right),j=-n+1,...,-\mathrm{1,0,1},...,n-1.$所得功率谱为互功率谱。本实施例的其他内容与实施例1相同。

上述计算完毕后，显示终端的显示如如图7和图8所示。