基于傅里叶变换的数字化脉冲对称高斯成形算法设计

引用本文

吴和喜, 陈玺婷, 袁新宇, 李思志, 刘义保. 基于傅里叶变换的数字化脉冲对称高斯成形算法设计[J]. 光谱学与光谱分析, 2023,43(8): 2483-2487.
WU He-xi, CHEN Xi-ting, YUAN Xin-yu, LI Si-zhi, LIU Yi-bao. Design for a Symmetrical Gaussian Shaping Algorithm Based on Fourier Transform for Handling Digital Pulse[J]. Spectroscopy and Spectral Analysis, 2023,43(8): 2483-2487.
Doi:10.3964/j.issn.1000-0593(2023)08-2483-05 复制到剪切板

Permissions

《光谱学与光谱分析》期刊社所有

基于傅里叶变换的数字化脉冲对称高斯成形算法设计

吴和喜¹, 陈玺婷¹, 袁新宇², 李思志¹, 刘义保¹

1.核技术应用教育部工程研究中心(东华理工大学), 江西南昌 330013

2.江西省辐射环境监督站, 江西南昌 330077

作者简介: 吴和喜, 1985年生, 核技术应用教育部工程研究中心(东华理工大学)教授 e-mail: wujing-7758@163.com

收稿日期: 2022-07-19 修回日期: 2022-10-10 接受日期: 2022-10-10

基金: 国家自然科学基金项目(12265003), 江西省自然科学基金项目(20212BAB201005) , 国家重点研发计划项目(2022YFC2807402)资助

摘要

随着高速模数转换器的出现, 数字化脉冲滤波成形算法被广泛应用于核仪器。现行的滤波成形算法中梯形成形算法输出脉冲宽度可调, 有利于提升核仪器的脉冲通过率; 高斯成形算法具有优越的信噪比, 使得核仪器的能量分辨率较优。为保留高斯成形算法的优点, 提升其对重叠核脉冲的分辨能力, 必须使其输出脉冲宽度变窄。该研究在假设输出脉冲形状为高斯函数、幅度与输入脉冲相等的情况下, 基于傅里叶变换推导了数字化脉冲对称高斯成形算法的递推公式。该递推公式仅由乘法和加法运算构成, 易于在现场可编程门阵列中实现。通过对理想脉冲仿真发现: 随着输出脉冲标准差的增加, 对称高斯成形算法成形后脉冲越宽, 对重叠脉冲的分辨能力也越弱; 对称高斯成形算法成形后脉冲宽度比梯形成形后最窄脉冲更小, 说明对称高斯成形算法具有更优的重叠脉冲分辨能力。通过对18种信噪比的脉冲各10 000次仿真发现: 在同一输出脉冲标准差下, 随着输入脉冲的信噪比增加, 输出脉冲幅度均值先增加, 当信噪比大于55 dB后输出脉冲幅度维持稳定; 随着输出脉冲标准差增加, 输出脉冲更平滑, 结合频谱图发现切断频率亦更小, 说明其去除噪声的能力更强, 有利于提升核仪器的能量分辨率。在全反射X荧光仪硬件平台上现场测试发现: 随着输出脉冲标准差增加, 谱仪能量分辨率提升了5.42%, 进一步证实了增大输出脉冲标准差可提升噪声去除能力、脉冲幅度提取亦更准确; 能谱总计数率下降了13.02%, 进一步证实了成形后脉冲越宽丢弃的重叠脉冲亦越多。上述结论充分证实本文所设计的对称高斯成形算法提升了高斯成形算法的重叠脉冲分辨能力, 并能有效应用于数字化多道脉冲幅度分析器实现脉冲高速实时处理。

关键词: 对称高斯成形; 重叠脉冲; 实时处理

中图分类号:O657.3 文献标志码:A

Design for a Symmetrical Gaussian Shaping Algorithm Based on Fourier Transform for Handling Digital Pulse

WU He-xi¹, CHEN Xi-ting¹, YUAN Xin-yu², LI Si-zhi¹, LIU Yi-bao¹

1. Engineering Research Center of Nuclear Technology Application (East China University of Technology), Ministry of Education, Nanchang 330013, China

2. Radiation Environment Supervision Station of Jiangxi Province, Nanchang 330077, China

Abstract

With the emergence of high-speed analog-to-digital converters, shaping algorithms for handling digital pulses was applied broadly in nuclear instruments. Trapezoidal shaper and Gaussian shaper are the common shaping algorithms. The former can increase the pulse passing rate due to its narrowed output pulse. The latter can make better energy resolution of nuclear instruments because of its attractive signal-to-noise ratio (SNR). The output gaussian shaper pulse has to narrow to add the discernment ability to piled-up pulses. The recurrence formula of the symmetrical gaussian shaping algorithm (SGSA) is deduced based on Fourier transform in the assumption that the shape of the output pulse obeys the gaussian function and its amplitude is equal to the corresponding input pulse. Moreover, its recurrence formula is composed of only multiplier and adder, which are easily implemented in field programmable gate arrays. After dealing with ideal pulses with SGSA, the result shows that the discernment ability to piled-up pulses is weaker along with the bigger σ (standard deviation) of its output pulse. Compare the shape of the output pulse of the trapezoidal shaper at the smallest rise and flat-top time with this method under σ=3 T_s. This proves this method has a better discernment ability to piled-up pulses. After dealing 180 000 simulating pulses under different SNR with SGSA, the result reflects that the average extractive amplitude increases at SNR<55 dB, and then maintain stability. Meanwhile, the output pulse is smoother after the bigger value of σ. Their frequency spectrums also prove it. The cutoff-frequency decreases, and the denoising effect rises when σ increases. It can increase the energy resolution of nuclear instruments. Applying the total-reflection X-ray fluorescence spectrometer shows that the energy resolution increases 5.42 % and total count rate decreases 13.02% when σ increases. The experiments further strengthen the above-simulating conclusion. Those results proved that SGSA improves the discernment ability to piled-up pulses and can be implemented high-speed real-time processing for a pulse on a multi-channel analyzer.

Keyword: Symmetrical gaussian shaping algorithm; Piled-up pulses; Real-time processing

文章图片

引言

1972 年, Radeka开始使用数字滤波成形算法处理核辐射探测器输出的数字化脉冲^[1], 但受限于当时的硬件水平, 难以实现实时处理。随着高速模数转换器(analog-to-digital converter, ADC)的出现, 数字化脉冲分析技术被广泛应用于核仪器。数字化脉冲滤波成形算法是数字化脉冲分析技术的重要组成部分。目前易于在现场可编程门阵列(field programmable gate array, FPGA)或数字信号处理器(digital signal processing, DSP)中实现数字化脉冲实时处理的滤波成形算法有高斯成形算法和梯形成形算法^[2]。传统的高斯成形算法多源于对模拟电路响应函数的离散化^{[3, 4]}, 基于电容积分特性将双指数脉冲变换成类高斯形状, 难以转变成真正的高斯函数形状并实现宽度可调。相比较而言, 高斯成形算法具有优越的信噪比^{[3, 5]}, 输出能谱的分辨率较优; 梯形成形算法输出的脉冲可调至更窄^[6], 对重叠核脉冲的分辨能力更强, 脉冲通过率更高。为兼顾高斯成形算法和梯形成形算法的优点, 本文拟基于傅里叶变换思想, 设计一款能输出满足高斯函数分布且宽度可调的数字化脉冲高斯成形算法(以下简称“ 对称高斯成形算法” ), 并通过仿真与现场测试验证所设计算法的可行性。

1 理论依据

对称高斯成形算法的原理是通过傅里叶变换将服从式(1)分布的输出脉冲与双指数分布的输入脉冲频谱之比变换至时域。

$y (t) = F e^{- t^{2} / (2 σ^{2})}$ (1)

式(1)中, y(t)为t时刻输出脉冲的电压值; F为输出脉冲的幅度值; σ 为高斯函数的标准差。

核辐射探测器收集到的电流信号经C-R微分电路后脉冲服从双指数分布^[2]

$x (t) = A (e^{- t / τ_{1}} - e^{- t / τ_{2}}) κ (t)$ (2)

式(2)中, x(t)为t时刻输入脉冲的电压值; κ (t)为阶梯函数, 在t≥ 0时κ (t)=1、其余时刻κ (t)=0; A、 τ ₁和τ ₂为待拟合参数。欲使输入脉冲x(t)与输出脉冲y(t)幅度相等, 则

$F = A (e^{- t_{F} / τ_{1}} - e^{- t_{F} / τ_{2}})$ (3)

式(3)中, 参数t_F值如式(4)

$t_{F} = τ_{1} τ_{2} (\ln τ_{2} - \ln τ_{1}) / (τ_{2} - τ_{1})$ (4)

根据傅里叶变换思想可知

$y (t) = \frac{1}{2 π} \int_{-}^{\infty} Y (ω) e^{jωt} d ω$ (5)

式(5)中, j为虚数基本单位, j²=-1; Y(ω )为输出脉冲y(t)的频谱。

假设输出脉冲频谱与输入脉冲频谱之比为H(ω ), 其表达式为

$H (ω) = Y (ω) / X (ω)$ (6)

式(6)中, X(ω )为输入脉冲x(t)的频谱。将式(6)代入式(5)可得

$y (t) = \frac{1}{2 π} \int_{-}^{\infty} X (ω) H (ω) e^{jωt} d ω$ (7)

假设h(T)为H(ω )在时域的函数, 即

$H (ω) = \int_{-}^{\infty} h (T) e^{- jωt} d T$ (8)

将其代入式(7)可得

$y (t) = \frac{1}{2 π} \int_{-}^{\infty} X (ω) e^{jωt} [\int_{-}^{\infty} h (T) e^{- jωT} d T] d ω$ (9)

改变积分先后顺序, 可得

$y (t) = \int_{-}^{\infty} h (T) [\frac{1}{2 π} \int_{-}^{\infty} X (ω) e^{- jω (t - T)} d ω] d T$ (10)

由x(t)与X(ω )间关系, 式(10)可化简为

$y (t) = \int_{-}^{\infty} h (T) x (t - T) d T$ (11)

可以看出, 将h(T)函数与输入脉冲x(t)卷积即可得到服从对称高斯分布的输出脉冲y(t)。根据傅里叶变换可知

$X (ω) = A \frac{τ_{1} - τ_{2}}{(1 + jω τ_{1}) (1 + jω τ_{2})}$ (12)

$Y (ω) = \sqrt[]{2 π} σF e^{- ω^{2} σ^{2} / 2}$ (13)

结合式(3)、式(6)、式(12)和式(13)可得

$H (ω) = \sqrt[]{2 π} σC [1 - ω^{2} τ_{1} τ_{2} + jω (τ_{1} + τ_{2})] e^{- ω^{2} σ^{2} / 2}$ (14)

式(14)中 $C = (e^{- t_{F} / τ_{1}} - e^{- t_{F} / τ_{2}}) / (τ_{1} - τ_{2}) 。由 h (T) 与 H (ω)$ 间关系可知

$h (T) = \frac{1}{2 π} \int_{-}^{\infty} H (ω) e^{jωT} d ω$ (15)

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

$h (T) = C \{1 - \frac{1}{σ^{2}} [T (τ_{1} + τ_{2}) + τ_{1} τ_{2} (1 - \frac{T^{2}}{σ^{2}})]\} e^{- T^{2} / (2 σ^{2})}$ (16)

因高速ADC以采样周期Ts将模拟脉冲离散成数字化脉冲, 为此在FPGA内运用式(11)处理前需先将其离散化, 具体如式(17)

$y (n) = \overset{N}{\sum_{i = - M}} h (i) x (n - i)$ (17)

式(17)中, i=T/T_s; n=t/T_s; x(n)和y(n)分别为t时刻高速ADC对x(t)和y(t)信号的采样值; h(i)为T时刻h(T)函数的值; M和N为正整数, 值的大小与计算精度、时长有关。

2 实验部分

2.1 理想脉冲仿真成形

在采样周期T_s=20 ns时对全反射X荧光仪中硅漂移探测器(silicon drift detector, SDD)输出脉冲进行数字化, 利用式(2)在拟合优度R²=0.982下拟合得到参数τ ₁=1 279.6 ns和τ ₂=40.2 ns, 代入式(4)可得t_F=143.6 ns。假设σ 值分别为3T_s、 5T_s和10T_s, 与τ ₁、 τ ₂和t_F值一并代入式(16)得到h(T)随T的变化规律如图1所示。并以采样周期T_s将h(T)离散化, 得到h(i)如图1中离散点所示。可以看出, 随着σ 值增大, h(T)函数向外延展的越宽, h(i)的取值也需越多, 算法的计算量亦越大。

	Figure Option View Download New Window
	图1 不同σ 值下h(T)函数分布Fig.1 The distribution of h(T) under different σ

将τ ₁和τ ₂以及A=1 785.6 mV代入式(2), 并在采样周期T_s下将x(t)数字化成x(n), 与图1所示的h(i)一起代入式(17)得到成形后的脉冲, 结果如图2所示。为减小计算量, 当|h(i)|< 0.01时令h(i)=0。可以看出: 所有σ 值下输出脉冲幅度与输入脉冲幅度相等且呈高斯函数分布; σ 值越小输出脉冲宽度Tσ (约为6σ /T_s)越窄, 且均小于输入脉冲宽度, 说明对称高斯成形算法具备较强的重叠脉冲分辨能力。与图2中梯形成形后最窄脉冲(上升沿与平顶持续时间均最短)比对, 发现σ =3T_s时对称高斯成形算法输出脉冲更窄, 说明对称高斯成形算法在重叠脉冲分辨能力上优于梯形成形算法。

	Figure Option View Download New Window
	图2 理想脉冲仿真成形结果Fig.2 Transform results of ideal pulse from simulation

2.2 含噪声脉冲仿真成形

因探测器输出脉冲上叠加了电子学噪声, 为评价对称高斯成形算法对含噪声信号的处理能力, 仿真了信噪比(signal-noise ratio, SNR)从15 dB到100 dB(间隔5 dB)共18种含噪声脉冲分别在σ =3T_s、 5T_s和10T_s下的成形结果, 典型情况如图3所示。可以看出: 成形后的输出脉冲呈现高斯函数分布; 随着σ 值增加成形后脉冲更为平滑, 说明其噪声过滤能力随σ 值增加而提升。将上述三种σ 值分别代入式(14), 得到对称高斯成形算法的频谱响应曲线如图4所示。可以看出: 对称高斯成形算法在频域展现为低通滤波器; 随着σ 值增加截断频率越小, 对称高斯成形算法对噪声的过滤能力亦越强; 与文献[4]中最佳参数时的数字化Sallen-Key滤波器相比, 对称高斯成形算法的截断频率更小噪声过滤能力更好。

	Figure Option View Download New Window
	图3 典型含噪声脉冲仿真成形结果Fig.3 Transform results of typical pulse including noise by simulation

	Figure Option View Download New Window
	图4 不同σ 值下的频率响应曲线Fig.4 The curves of frequency response under different σ

将上述18种含噪声脉冲各仿真10 000次, 在σ =5T_s下分别成形, 提取最大值点作为脉冲幅度值, 典型(SNR=55 dB)含针声脉冲幅度提取结果如图5所示。可以看出10 000次脉冲幅度提取值服从正态分布, 在拟合优度R²=0.993下得到其平均值为2 168.74 mV、标准差为1.34。据此计算出的上述18种含噪声脉冲仿真成形后幅度提取值的变化范围如图6所示, 可以看出: 随着SNR值增加, 幅度均值先上升后维持稳定(SNR> 55dB); 叠加的噪声越弱, 成形后幅度变化范围降低, 对谱仪的能量分辨率影响亦减小。

	Figure Option View Download New Window
	图5 典型含噪声脉冲幅度提取结果Fig.5 The extractive amplitudes of typical pulse including noise

	Figure Option View Download New Window
	图6 含不同SNR噪声脉冲成形后幅度提取值Fig.6 The extractive amplitudes of pulse under different SNR

3 结果与讨论

实验用全反射X荧光仪硬件结构如文献[7]所示, SDD输出脉冲进入以AD8066集成放大器为核心的信号调理电路, 实现限幅处理与信号放大。因ADC位数越高, 量化误差越小^[6]。采用14位高速ADC芯片将信号调理电路输出模拟脉冲于采样周期T_s=20 ns转换成数字化脉冲, 并存入FPGA内核心芯片EPCS16的闪存FIFO中(每次存入连续1 024个ADC采样值), 典型数字化脉冲如图7(a)所示, 在FPGA内按照图8所示流程加以处理, 具体如下:

	Figure Option View Download New Window
	图7 实测连续脉冲成形结果Fig.7 Transform results of a series of digital pulse with different shaping Algorithm

	Figure Option View Download New Window
	图8 FPGA内脉冲处理流程Fig.8 Process for handling pulse in FPGA

(1) 基线计算: 由于存储在FIFO中的x(n)数据序列时长仅为20.48 μ s, 如此短时间内的基线变化可忽略。基线值B计算采用最小平均值法^[8], 即提取相对最小的8个x(n)值求平均。

(2)脉冲识别: 对x(n)进行冲击成形(算法详见文献[7]), 结果如图7(b)所示, 可以看出各脉冲起始位置出现一个幅度与之相等的冲击响应, 而其他位置变换值均相对较小, 通过阈值(本文阈值取10)甄别即可确定任意第j个脉冲的起始位置S_j。

(3)幅度提取: 在给定σ 值下利用式(17)对x(n)进行变换得到y(n), 结果如图7(c)— (e)所示, 可以看出: 当σ =10T_s时成形结果无法有效区分前6个重叠脉冲; 当σ =5T_s时可有效区分出第4~10个重叠脉冲, 但依然无法区分前3个重叠脉冲; 当σ =3T_s时可有效区分出所有的10个重叠脉冲。因为随着σ 值减小成形后核脉冲宽度越窄, 所以能够区分的重叠脉冲间隔亦越窄。以第j个脉冲为例, 如果S_j₊₁-S_j< T_σ, 说明成形后脉冲宽度小于相邻脉冲间隔, 可有效提取第j个脉冲幅度值, 即F_j=max[y(n), S_j< n< S_j₊₁]; 反之则丢弃第j个脉冲不记录。

(4)基线恢复: 将提取的第j个脉冲幅度减去基线值并向下取整, 得到第j个脉冲的实际幅度值floor(F_j-B)。

(5)能谱记录: 将第floor(F_j-B)道存储地址内数据加1。

依据步骤(1)— 步骤(5)记录能有效提取幅度的脉冲。继续获取1 024个ADC采样值以上述过程一一分析, 直至测量时间结束。

分别在σ =3T_s、 5T_s和10T_s下多次测量Fe的能谱, 应用GammaVision软件计算能谱总计数率与能量分辨率, 结果如表1所示。从表1可以看出: 随着σ 值增加, 谱仪能量分辨率提升了5.42%, 说明去除了的噪声越多脉冲幅度提取准确度越高; 能谱总计数率下降了13.02%, 说明成形后脉冲越宽丢弃的重叠脉冲越多。

表1 不同σ 值下铁样品多次实测谱分析结果 Table 1 Analyzing results from many measuring spectrums of iron sample under different σ

4 结论

在已知时变系统输入与输出脉冲形状的基础上, 基于傅里叶变换法推导了数字化脉冲对称高斯成形算法的递推公式。对特定谱仪而言, 对称高斯成形算法仅与输出脉冲半高宽有关。通过仿真发现随着输出脉冲半高宽增加, 对称高斯成形算法的平滑效果越好, 但对重叠脉冲分辨能力会变差。将其应用于全反射X荧光仪现场实验测试后发现: 随着设置的输出脉冲半高宽增加, 能谱仪的总计数率下降了13.02%、能量分辨率提升了5.42%。后续应用中应根据测量要求调整输出脉冲半高宽, 并将对称高斯成形算法应用于其他核辐射探测器中充分检验其普适性。

参考文献

文献列表

[1]	Radeka V. Nuclear Instruments and Methods, 1972, 99(3): 525. [本文引用:1]
[2]	Glenn F. Knoll. Radiation Detection and Measurement. 4th ed. New York: John Wiley & Sons, Inc, 2010. [本文引用:2]
[3]	Liu Yinyu, Zhang Jinglong, Liu Lifang, et al. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 2018, 906: 1. [本文引用:2]
[4]	Zhang Huaiqiang, Tang Bin, Wu Hexi, et al. Nuclear Science and Techniques, 2019, 30(10): 145. [本文引用:2]
[5]	Du Siyan, Liu Jia, Yang Fengfan, et al. Energy Reports, 2020, 6(Supplement 9): 343. [本文引用:1]
[6]	SONG Wen-gang, ZHANG Li-jun, ZHANG Jing, et al(宋文刚, 张立军, 张晶, 等). Acta Physica Sinica(物理学报), 2022, 71(1): 012903. [本文引用:2]
[7]	WU He-xi, DI Run-jie, LIU Yu-juan, et al(吴和喜, 邸润洁, 刘玉娟, 等). Spectroscopy and Spectral Analysis(光谱学与光谱分析), 2021, 41(7): 2148. [本文引用:2]
[8]	CHEN Wei, ZHOU Jian-bin, FANG Fang, et al(陈伟, 周建斌, 方方, 等). Nuclear Techniques(核技术), 2018, 41(5): 050401. [本文引用:1]

1972

0.0

... 引言1972 年, Radeka开始使用数字滤波成形算法处理核辐射探测器输出的数字化脉冲^[1], 但受限于当时的硬件水平, 难以实现实时处理 ...

2010

0.0

... 目前易于在现场可编程门阵列(field programmable gate array, FPGA)或数字信号处理器(digital signal processing, DSP)中实现数字化脉冲实时处理的滤波成形算法有高斯成形算法和梯形成形算法^[2] ...

... 核辐射探测器收集到的电流信号经C-R微分电路后脉冲服从双指数分布^[2] ...

2018

0.0

... 传统的高斯成形算法多源于对模拟电路响应函数的离散化^[3,4], 基于电容积分特性将双指数脉冲变换成类高斯形状, 难以转变成真正的高斯函数形状并实现宽度可调 ...

... 相比较而言, 高斯成形算法具有优越的信噪比^{[3, 5]}, 输出能谱的分辨率较优 ...

2019

0.0

... 与文献[4]中最佳参数时的数字化Sallen-Key滤波器相比, 对称高斯成形算法的截断频率更小噪声过滤能力更好 ...

2020

0.0

... 相比较而言, 高斯成形算法具有优越的信噪比^{[3, 5]}, 输出能谱的分辨率较优 ...

2022

0.0

... 梯形成形算法输出的脉冲可调至更窄^[6], 对重叠核脉冲的分辨能力更强, 脉冲通过率更高 ...

... 因ADC位数越高, 量化误差越小^[6] ...

2021

0.0

... 3 结果与讨论实验用全反射X荧光仪硬件结构如文献[7]所示, SDD输出脉冲进入以AD8066集成放大器为核心的信号调理电路, 实现限幅处理与信号放大 ...

... (2)脉冲识别: 对x(n)进行冲击成形(算法详见文献[7]), 结果如图7(b)所示, 可以看出各脉冲起始位置出现一个幅度与之相等的冲击响应, 而其他位置变换值均相对较小, 通过阈值(本文阈值取10)甄别即可确定任意第j个脉冲的起始位置S_j ...

2018

0.0

... 基线值B计算采用最小平均值法^[8], 即提取相对最小的8个x(n)值求平均 ...