-
高速采样率下单位时间内数据量巨大,因为无法实时处理会舍弃部分数据,导致容易漏掉偶发的瞬态信号[1]。为了提高瞬态信号的捕获概率,需要借助瞬态信号检测技术,定位和捕获感兴趣的信号成分。传统信号检测技术利用信号时域信息来检测,如信号的幅度、上升沿或下降沿、脉冲宽度等[2–4]。然而,传统时域检测方法无法从信号频域对信号进行检测。如现代雷达系统测试需要监测频谱的分配情况[5]、现代频谱管理需要监测选定频带中的瞬态干扰或脉冲噪声[6]。可见,研究以目标频域为对象的检测方法具有重要意义。通过频域特征来检测信号成为了测量领域的重要研究方向,被广泛应用在测试仪器的频域触发系统中,实现的触发方式包括频率边沿触发、频率模板触发(FMT)和DPX密度触发[7–10]。
然而,这些触发方法对目标频段瞬态信号的检测尚存在局限。频率边沿触发利用相位变化率(瞬时频率)来定位信号的频率变化。但它是时间的单值函数,在给定时间内只能给出一个瞬时频率值,当输入是多音信号时,信号的相位信息是多个频率分量综合作用的结果,不能从多分量信号中分离出感兴趣的频率成分。为了分离不同的频率成分,可采用频率模板触发和DPX密度触发。频率模板触发根据用户定义的目标频段和目标振幅来查找和定位用户感兴趣的波形,DPX密度触发机制是根据目标频段和目标振幅,检测超过统计阈值的信号事件。但当信号是持续时间较短的瞬态信号时,DFT无法提供足够的频率分辨率来感知信号能量分布,从而导致频谱泄漏[11]。此外,对于经常遇到的分析窄带的应用场景,带外频谱计算是多余的。而使用DFT计算整个奈奎斯特域的频谱,冗余的计算会降低信号检测速度。因此,现有的频域信号检测方法无法可靠地解决目标频段瞬态信号检测问题,有必要提出一种更高效率和更可靠的检测方法来捕获目标频段中的瞬变事件。
本文提出了一种基于Slepian 序列信号字典的检测方法。首先选取标准正交的Slepian序列来组成信号字典,其包含了目标频段内信号的所有波形特征,其次通过判断被观测信号与字典的匹配程度,来定位目标频段内的瞬态信号。
-
被分析的信号
$x(t)$ 是时间的连续函数。数字系统对模拟信号进行离散化处理,得到信号样点。由于数字系统中乘法器、加法器等硬件资源限制,单次只能分析处理有限的数据量,如图1所示。该数据量记为${N_F}$ ,则信号的观测时间2T为:$$ 2T = {N_F}{T_s} $$ (1) 式中,
${T_s}$ 表示离散抽样间隔。相应的采样数据可表示成$ {x_n} = {{x(}}n{T_s}{\text{)}},n = 0,1, \cdots ,{N_F} - 1 $ 。假设目标频带是基带
$ [ - \varOmega ,\varOmega ] $ ,$\varOmega $ 表示角频率(单位:rad/s),与信号频率$f$ (单位:Hz)和采样率${f_s}$ (单位:Hz)有关:$$ \varOmega {\text{ = 2}}{\text{π}} f/{f_s} $$ (2) 该频带内信号能量可由采样数据
$ {x_n} $ 估计:$$ \int\nolimits_{ - \varOmega }^\varOmega {{{\left| {X\left( \xi \right)} \right|}^2}\frac{{{\mathrm{d}}\xi }}{{2{\text{π}} }}} = \sum\limits_{n = 0}^{{N_F} - 1} {\sum\limits_{m = 0}^{{N_F} - 1} {{x_n}{x_m}\frac{{\sin [\varOmega (n - m)]}}{{{\text{π}} (n - m)}}} } $$ (3) 式中,
$X(\xi )$ 表示$\{ {x_n}\} $ 的傅里叶变换:$$ X\left( \xi \right) = \sum\limits_{n = 0}^{{N_F} - 1} {{x_n}{{\mathrm{e}}^{ - {\mathrm{j}}\xi n}}} \quad\xi \in [ - {\text{π}} ,{\text{π}} ] $$ (4) 将
$\{ {x_n}\} $ 看作向量,则其所在的向量空间维数为${N_F}$ 。Laudau、Pollak和Slepian找到了该空间的一组标准正交基,称作Slepian序列$ \{ \upsilon _n^{\left( k \right)},1 \leqslant k \leqslant {N_F}\} $ [12–16]:$$ \left\langle {\upsilon _n^{\left( {{k_1}} \right)},\upsilon _n^{\left( {{k_2}} \right)}} \right\rangle = \sum\limits_{n = 0}^{{N_F} - 1} {\upsilon _n^{\left( {{k_1}} \right)}\overline {\upsilon _n^{\left( {{k_2}} \right)}} } = \left\{ {\begin{array}{*{20}{c}} 1&{{k_1} = {k_2}} \\ 0&{{k_1} \ne {k_2}{\text{}}} \end{array}} \right. $$ (5) 用
$ \upsilon _n^{\left( k \right)} $ 替换式中的${x_n}$ ,同样可计算Slepian序列位于目标频带的能量${\lambda _k}$ :$$ {\lambda _k} = \sum\limits_{n = 0}^{{N_F} - 1} {\sum\limits_{m = 0}^{{N_F} - 1} {\upsilon _n^{\left( k \right)}\upsilon _m^{\left( k \right)}\frac{{\sin [\varOmega (n - m)]}}{{{\text{π}} (n - m)}}} } $$ (6) ${\lambda _k}$ 取值只取决于乘积$T\varOmega $ ,并随$k$ 增大而减小,如图2所示。阶数$k$ 靠前的序列其能量与1接近,当$k$ 超过阈值${\text{2}}T\varOmega /{\text{π}} $ ,能量迅速衰减到0,之后停留在0附近。衰减区域的宽度与$\log (T\varOmega )$ 成正比。基序列在目标频带被观测的能量${\lambda _k}$ 越大,说明其特征与目标频带的信号特征越相关。能量衰减大的序列与目标频带的信号特征不符,对检测毫无用处。图3给出了在$N_{\mathrm{F}}=83,\;2T\varOmega/\text{π}=2.5 $ 情况下,${\lambda _k}$ 接近1的Slepian序列波形。 -
将与目标频带的信号特征相关的所有Slepian序列集合起来,构成信号字典
$ D = \{ \upsilon _n^{\left( k \right)},1 \leqslant k \leqslant K\} $ ,这里$K = 2T\Omega /{\text{π}} $ 。由字典内的序列张成的空间为$V$ ,具备有限的时间和频率支撑,位于图4中的$ [ - T,T] \times [ - \varOmega ,\varOmega ] $ 区域。$V$ 的自由度对应字典里的序列数量。本文方法通过判断被测信号与字典的时频特征是否一致,实现对信号的检测。总体检测方案如图5所示,检测系统根据设定的
${N_F}$ 和目标带宽$\Delta f$ 确定字典包含的序列个数K,结合式(1)和式(2),有:$$ K = 2T\varOmega /{\text{π}} = {N_F}\Delta f/{f_s} $$ (7) 根据
${N_F}$ 和$\Delta f$ 生成相应的Slepian序列,生成方法参见文献[14]。注意到
$V$ 的时频特征位于基带。当目标频带区域位于图4中的中心频率${\varOmega _c}$ 时,还需要时间移位和频率移位,即图5中对信号处理部分,采集的数据样点存储在图5中的环形数据缓冲器内。时间移位通过不同的起始读地址实现,将时间中心位置不同的数据帧读出,读地址取决于当前的计数器输出的起始样点位置,起始点计数器的计数步进Z可调。读出的数据帧$ \{ {x_n}\} $ 以步进Z滑动,如图6所示。帧1样点与信号无关,帧2有部分样点来自信号,帧3样点全部来自信号。因此,帧3定位到了瞬态信号。为了提高时间定位的准确性,帧步进$Z$ 应设置较小,这时数据帧之间存在重叠。但$Z$ 越小,检测的计算复杂度越大。频率移位可通过数字下变频实现,设定可调本振的频率与目标带宽的中心频率${\varOmega _c}$ 一致。其具体实现如图7所示,输入数据${x_n}$ 与两路正交本振混频,输出同相分量${I_n} = {x_n} \cos ({\varOmega _c}n)$ 和正交分量${Q_n} = - {x_n}\sin ({\varOmega _c}n)$ 。两路正交本振信号由数控振荡器NCO和相应的移相器产生。如图4,经过数据预处理,目标时频特征与信号字典的特征重合,即预处理后的信号
${u_n}$ 为:$$ {u_n} = {I_n} + {\mathrm{j}}{Q_n} $$ (8) 借助正交分解,当
${u_n}$ 的特征与字典匹配时,${u_n}$ 在空间$V$ 中的投影能量达到最大。所以,投影能量衡量着${u_n}$ 与字典的匹配度。为了便于计算,将${u_n}$ 由字典中的向量$\boldsymbol \upsilon _n^{\left( k \right)}$ 线性表出:$$ {u_n} = \sum\limits_{k = 1}^{{\text{2}}T\varOmega /{\text{π}} } {\left\langle {{u_n},\boldsymbol \upsilon _n^{\left( k \right)}} \right\rangle } \boldsymbol \upsilon _n^{\left( k \right)} + {\mathrm{R}}{{\mathrm{u}}_n} $$ (9) 式中,
$ {\mathrm{R}}{{\mathrm{u}}_n} $ 表示${u_n}$ 结构与字典向量不匹配的残留项(Residue)。${u_n}$ 在函数空间$V$ 中投影(Projection)的能量表示为:$$ {\left\| {{P_V}{u_n}} \right\|^2} = \sum\limits_{k = 1}^{{\text{2}}T\varOmega /{\text{π}} } {{{\left| {\left\langle {{u_n},\boldsymbol \upsilon _n^{\left( k \right)}} \right\rangle } \right|}^2}} $$ (10) 用功率
${P_L}$ (dBm)表示的匹配度为:$$ {P_L} = 10{\log _{10}}\left( {\frac{{{{{{\left\| {{P_V}{u_n}} \right\|}^2}} \mathord{\left/ {\vphantom {{{{\left\| {{P_V}{u_n}} \right\|}^2}} {{N_F}}}} \right. } {{N_F}}}}}{{{Z_L}}}} \right) $$ (11) 式中,
${Z_L}$ 表示50 Ω载阻抗。图5中,将计算出的匹配度${P_L}$ 与比较器相接,当大于设定的阈值,输出触发信号,指示目标信号的出现。 -
为了说明本文方法的高效性,与短时傅里叶变换(STFT)、离散小波变换(DWT)和加窗Wigner-Ville分布(Pseudo WVD, PWVD)进行对比。对比实验假定采样率为122.88 MHz,采样数据总长度为
$N$ ,系统单次处理数据帧长度为${N_F}$ 。表1为不同分析方法的参数设置。表 1 对比实验中的参数设置
对于所提出的方法,构建信号字典和数据预处理,可在分析前完成,在分析过程中不产生计算操作。计算复杂度来自于计算匹配度
$ {\left\| {{P_V}{x_n}} \right\|^2} $ ,对于单次数据帧分析,需要K次乘累加运算,计算${u_n}$ 在$\boldsymbol \upsilon _n^{\left( k \right)}$ 投影。乘累加运算量${N_F}$ 为:$$ \lt {u_n},\boldsymbol \upsilon _n^{(k)} \gt = \sum\limits_{n = 0}^{{N_F} - 1} {{u_n}\boldsymbol \upsilon _n^{(k)}} \quad k = 1,2, \cdots ,K $$ (12) 需要K次加法将K个投影的能量求和。假设截取窗口滑动步进为
$Z$ ,则处理的数据帧总量$M = {\mathrm{floor}}\left( {\left( {N - {N_F}} \right)/Z} \right) + 1$ ,上述的计算要重复M次完成整个采样数据的分析。由此,提出的检测方法所需计算复杂度为:$$ {C_{{\text{Slepian}}}} = M\left( {K{N_F} + K} \right) $$ (13) 对于STFT,可使用快速傅里叶变换(Cooley Tukey FFT)算法,节省加法和乘法次数。假设窗函数
$w\left( k \right)$ 长度为${N_F}$ ,滑动步进$Z$ ,那么要处理的数据帧总量与本文提出的方法相同,其计算复杂度为:$$ {C_{{\text{STFT}}}} = M\left( {{N_F} + {N_F}{{\log }_2}{N_F}} \right) $$ (14) 对于DWT,也可借助FFT算法 [20]:
$$ {W_x} {{a_j},{b_m}} = {F_n}^{ - 1} {\left( {{F_n}x} \right) \odot {{{\widehat {\boldsymbol{\psi}} }_j}} } \quad j = 0, 1,\cdots ,J - 1 $$ (15) 式中,
${a_j} = {r^j}\Delta u$ 和${b_m} = m\Delta u$ 是对尺度和时移的离散化;$J$ 表示允许的频带数量;$ \odot $ 表示向量中元素相乘;向量$ {{\boldsymbol{\widehat \psi}} _j} $ 的长度与$N$ 一致,$ {({{\boldsymbol{\widehat \psi}} _j})_m} = {a_j}^{1/2}{\boldsymbol{\widehat \psi}} ({a_j}{\xi _m}) $ ,${\xi _m}$ 表示频率区间抽样($ {\xi _m} = 2{\text{π}} m/N, \; m = 0,1, \cdots ,N - 1 $ )。DWT的计算复杂度为:$$ {C_{{\text{DWT}}}} = J\left( {N{{\log }_2}N + N} \right) + N{\log _2}N $$ (16) 比较分析使用Morse小波作为
$ {{\boldsymbol{\widehat \psi}} _j} $ ,该类小波适宜瞬态信号分析,$J$ 取59[18]。对于PWVD,其计算复杂度为[19]:
$$ {C_{{\text{PWVD}}}} = {N^2}/2 + M\left( {2{N_w} + 2{N_w}{{\log }_2}2{N_w}} \right) $$ (17) 式中,
$2{N_w}$ 表示窗函数长度。图8给出在不同N情形下,采样率为12288 MHz时不同分析方法的计算复杂度。图中横坐标表示数据容量,纵坐标表示计算复杂度,以对数形式表示。随着
$N$ 的增加,各类方法的计算复杂度增大。在$256 \leqslant N \leqslant 4\;096$ 区间内,本文的方法比其他方法的计算复杂度低,计算复杂度最大的是STFT,其次是PWVD和DWT,PWVD的计算复杂度增长率快于DWT,在$N{\text{ = 1\;024}}$ 处,PWVD的计算复杂度超过DWT。图8还显示了本文方法在情形1的计算复杂度始终高于情形2。这是由于在
${N_F}W$ 和滑动步进$Z$ 相同的条件下,情形2的目标带宽$\Delta f$ 更宽,需处理的帧长度${N_F}$ 从61减少到了27,计算复杂度降低。定义计算量减少的百分比
$ \% {\mathrm{CR}} $ 为:$$ \text{%} {\mathrm{CR}} = 100 \text{%} \times \frac{{{C_X} - {C_{{\text{Slepian}}}}}}{{{C_X}}} $$ (18) 式中,
${C_X}$ 代表计算复杂度(${C_{{\text{STFT}}}}$ ,${C_{{\text{DWT}}}}$ 和${C_{{\text{PWVD}}}}$ )。根据定义,表2给出了$ \text{%} {\mathrm{CR}} $ 结果,说明本文方法的计算复杂度比STFT减少92%以上,比DWT减少71%以上,比PWVD减少35%以上。此方法是对传统时频分析方法的改良,即它只感知局部频谱,保留感兴趣波段的分析计算,从而使得计算复杂度降低,提高了信号检测效率。表 2 计算量减少百分比
% 方法 %CR,情形1/情形2 N=256 N=512 N=1024 N=2048 N=4096 STFT[17] 93.2/94.4 92.8/94.2 92.6/94.0 92.5/94.0 92.5/94.0 DWT[18] 71.7/76.7 72.7/77.8 74.3/79.3 76.1/80.8 77.7/82.1 PWD[19] 35.4/46.7 56.1/64.4 73.1/78.3 84.9/87.9 92.0/93.6
Transient Signal Detection Method for Target Frequency Band Using Slepian Series
-
摘要: 提出一种基于Slepian序列信号字典的目标频段瞬态信号检测方法,该方法只关注目标频段的能量信息,是一种高效率的检测方式。首先,选取标准正交的Slepian序列组成信号字典,该字典能够表征目标频段内信号特征;然后,通过判断被观测信号样点与字典的匹配程度实现检测。对比实验表明,该方法的计算效率比短时傅里叶变换提升92%以上,比离散小波变换提升71%以上,比加窗Wigner-Ville分布提升35%以上。仿真实验使用时域稀疏的脉冲调制信号进行验证,结果表明了该检测方法的有效性。Abstract: A huge amount of data will be generated per unit of time at high sampling rates. Some data will be discarded when the data cannot be processed in real-time, thus missing the episodic transient signals. Practical test scenarios are often faced with sparse signals in the frequency domain, and the concern is often narrow frequency bands. Therefore, a transient signal detection method for the target band is proposed based on a signal dictionary of the Slepian series, which is efficient by focusing only on target band energy information. First, a set of orthonormal Slepian sequences are selected to form a signal dictionary, which can characterize the signal features in the target band; then, the detection is achieved by judging the matching degree between the observed signal samples and the dictionary. The comparison experiment shows that the computational complexity of the proposed detection method is reduced by more than 92% compared with the short-time Fourier transform, more than 71% compared with the discrete wavelet transform, and more than 35% compared with the windowed Wigner-Ville distribution. Simulation experiments are performed using pulse-modulated signals with a sparse time domain for verification, and the results show the effectiveness of the proposed detection method. Based on the detection results, the sampling data that are not of interest can be discarded, which saves storage resources and reduces the amount of data for signal processing.
-
表 1 对比实验中的参数设置
表 2 计算量减少百分比
% 方法 %CR,情形1/情形2 N=256 N=512 N=1024 N=2048 N=4096 STFT[17] 93.2/94.4 92.8/94.2 92.6/94.0 92.5/94.0 92.5/94.0 DWT[18] 71.7/76.7 72.7/77.8 74.3/79.3 76.1/80.8 77.7/82.1 PWD[19] 35.4/46.7 56.1/64.4 73.1/78.3 84.9/87.9 92.0/93.6 -
[1] 刘洪庆, 向前. 示波器最新技术进展与发展趋势[J]. 电子质量, 2021(8): 1-5. doi: 10.3969/j.issn.1003-0107.2021.08.001 LIU H Q, XIANG Q. Latest technology progresses and development trends of oscilloscope[J]. Electronics Quality, 2021(8): 1-5. doi: 10.3969/j.issn.1003-0107.2021.08.001 [2] TEKTRONIX. Fundamentals of real-time spectrum analysis[EB/OL]. (2019-05-23). https://download.tek.com/document/37W_17249_6_Fundamentals_of_Real-Time_Spectrum_Analysis1.pdf. [3] TEKTRONIX. Fundamentals of the MDO4000 series mixed domain oscilloscope[EB/OL]. (2019-05-23). https://www.tek.com/en/documents/application-note/fundamentals-mdo4000-series-mixed-domain-oscilloscope. [4] 李毅, 师奕兵, 王厚军, 等. 数字存储示波器触发电路的数字化技术研究[J]. 仪器仪表学报, 2004, 25(3): 385-387. doi: 10.3321/j.issn:0254-3087.2004.03.027 LI Y, SHI Y B, WANG H J, et al. Research on the digitization technique of trigger circuit in digital storage oscilloscope[J]. Chinese Journal of Scientific Instrument, 2004, 25(3): 385-387. doi: 10.3321/j.issn:0254-3087.2004.03.027 [5] TEKTRONIX. Advanced radar analysis tools for measuring modern radars[EB/OL]. (2020-11-10). https://www.tek.com/en/documents/application-note/advanced-radar-analysis-tools-measuring-modern-radar-application-note. [6] TEKTRONIX. Spectrum management/surveillance solutions fact sheet[EB/OL]. (2019-05-23). https://www.tek.com/en/documents/fact-sheet/spectrum-management-surveillance-solutions-fact-sheet. [7] HILLMAN JR A K, HARWOOD S L. Time qualified frequency mask trigger: US8255179[P]. 2012-08-28. [8] STANTON S W, GEE E C, HILLMAN A K, et al. Frequency mask trigger with non-uniform bandwidth segments: US9702907[P]. 2017-07-11. [9] LIU B Q, WANG Q, YAN X, et al. A new spectrum measurement trigging method and its implementation[C]//Proceedings of the International Conference on Computational Problem-Solving. New York: IEEE, 2013, DOI: 10.1109/ICCPS.2013.6893538. [10] LI T, CHEN Z P. A wideband reconnaissance receiver design based on real-time spectrum analysis technology[C]//Proceedings of the IEEE 2nd Advanced Information Technology, Electronic and Automation Control Conference. New York: IEEE, 2017: 1463-1467. [11] 程佩青. 数字信号处理教程: MATLAB版[M]. 5版. 北京: 清华大学出版社, 2017. CHENG P Q. Digital signal processing course: MATLAB version[M]. 5th ed. Beijing: Tsinghua University Press, 2017. [12] KARNIK S, ROMBERG J, DAVENPORT M A. Thomson’s multitaper method revisited[J]. IEEE Transactions on Information Theory, 2022, 68(7): 4864-4891. doi: 10.1109/TIT.2022.3151415 [13] DAUBECHIES I. Ten lectures on wavelets[M]. Philadelphia, PA: Society for Industrial and Applied Mathematics, 1992. [14] BELL B, PERCIVAL D B, WALDEN A T. Calculating Thomson’s spectral multitapers by inverse iteration[J]. Journal of Computational and Graphical Statistics, 1993, 2(1): 119-130. doi: 10.1080/10618600.1993.10474602 [15] MOORE I C, CADA M. Prolate spheroidal wave functions, anintroduction to the Slepian series and its properties[J]. Applied and Computational Harmonic Analysis, 2004, 16(3): 208-230. doi: 10.1016/j.acha.2004.03.004 [16] HAYKIN S, THOMSON D J, REED J H. Spectrum sensing for cognitive radio[J]. Ingeniare Revista Chilena de Ingeniería, 2012, 20(2): 197-210. [17] 胡广书. 现代信号处理教程[M]. 北京: 清华大学出版社, 2004. HU G S. Modern signal processing course[M]. Beijing: Tsinghua University Press, 2004. [18] LILLY J M. Element analysis: A wavelet-based method for analysing time-localized events in noisy time series[J]. Proc Math Phys Eng Sci, 2017, 473(2200): 20160776. [19] GUNER K K, GULUM T O, ERKMEN B. FPGA-based wigner-hough transform system for detection and parameter extraction of LPI radar LFMCW signals[J]. IEEE Transactions on Instrumentation Measurement, 2021, 70: 3060584. [20] THAKUR G, BREVDO E, FUČKAR N S, et al. The Synchrosqueezing algorithm for time-varying spectral analysis: Robustness properties and new paleoclimate applications[J]. Signal Processing, 2013, 93(5): 1079-1094. doi: 10.1016/j.sigpro.2012.11.029 [21] JIA X, TANG L H, LIU S G, et al. Frequency-tunable pulsed microwave waveform generation based on unbalanced single-arm interferometer excited by near-infrared femtosecond laser[J]. Applied Sciences, 2021, 11(24): 11928. doi: 10.3390/app112411928