一种基于离散傅立叶变换的低成本窄带干扰频率估计方法与流程

1.本发明涉及一种窄带干扰频率估计方法，特别是一种基于离散傅立叶变换的低成本窄带干扰频率估计方法。


背景技术：

2.在通信技术和许多使用场景下，窄带干扰为常见的噪声干扰。设计上一般会使用陷波滤波器进行抑制或消除，而在作动此一滤波器前，会先判断该干扰是否为对讯号有严重影响，评量方式一般都以该干扰的能量是否大于讯号，或一个特定默认的容忍值，再针对估算后的干扰频率位置进行抑制或消除。
3.因一个陷波滤波器一次只能处理一个单频频率，芯片通常在成本和设计考虑下，所使用的陷波处理器数量会有其上限，假若陷波滤波器在处理nbi(narrow band interference，nbi，窄带干扰)时，消除干扰的频率位置有所偏差，会导致干扰消除不确实，使得部分干扰能量依然存在，最终仍将影响到通信效能。如果不能精确地找出干扰频率位置，使陷波滤波器能完整发挥其功用，将导致陷波滤波器无法效消除干扰源。
4.顾名思义，窄带干扰是一种占据相对窄的带宽的一种干扰。干扰的带宽(bw，bandwidth)与通信技术的所需信号的带宽进行比较。一个宽松的定义为，当干扰bw
interference
/bw
sigal
《0.01且功率谱密度(psd，power spectral density，功率谱密度)大于信号时，可以归类为窄带干扰nbi。窄带干扰可以进一步分为同频干扰和带外干扰，这取决于窄带干扰的频谱是否与通信信号重迭。这种窄带干扰对大部分的通信系统会有影响，严重时甚至于无法正常通信。在这种情况下，窄带干扰抑制或消除至关重要。
5.如图2所示，为窄带干扰与通信信号示意图，其中nbis为窄带干扰，signals为通信信号。
6.现行技术中已存在了许多nbi消除技术，包括预测技术、最小均方误差(mmse，minimum mean square error)检测器、变换域技术和陷波滤波。
7.在众多通信系统中，plc面对着许多的潜在干扰信号，尤其是在数百千赫兹以下的频率范围内，nbi的问题使得plc通信变得非常具有挑战性。根本上nbi不易排除可能的原因包括了：
8.1、plc通信系统频率范围内这些非线性负载的最大发射电平没有足够的标准。
9.2、开关式电源(来自plc模块的板载电源或电气设备/设备)在开关频率及其谐波频率上产生窄带干扰。
10.3、nbi功率通常很强，会压倒潜在的通信信号，因此难以建立可靠的通信链路。
11.4、nbi的存在是时变的。


技术实现要素：

12.发明目的：本发明所要解决的技术问题是针对现有技术的不足，提供一种基于离散傅立叶变换的低成本窄带干扰频率估计方法。
13.为了解决上述技术问题，本发明公开了一种基于离散傅立叶变换的低成本窄带干扰频率估计方法，包括以下步骤：
14.步骤1，对于输入信号进行离散傅里叶变换，得到输出频点；
15.步骤2，在所有输出频点中找到输出最大功率的频点作为频点索引；
16.步骤3，根据所述频点索引计算得到窄带干扰频率，完成基于离散傅立叶变换的低成本窄带干扰频率估计方法。
17.步骤2中所述的在所有输出频点中找到输出最大功率的频点作为频点索引k
peak
的方法包括：
[0018][0019]
其中，p(k)表示离散傅立叶变换dft输出频点k的频点功率值，d表示频点k的查找范围，n是离散傅立叶变换dft的大小。
[0020]
步骤3中所述的根据所述频点索引计算得到窄带干扰频率的方法包括：根据所述频点索引进行迭代搜索，具体方法如下：
[0021]
步骤3-1，定义降采样倍数l，确定搜索范围w；
[0022]
步骤3-2，在每一次迭代搜索过程中，将所述输入信号即离散时域信号样本直接降采样l倍，得到l倍的窄带干扰频率；
[0023]
步骤3-3，根据降采样后的最大功率的新频点，更新窄带干扰频率的估计值；
[0024]
步骤3-4，重复步骤3-2至3-3的迭代搜索过程，直至达到迭代终止条件，完成迭代搜索。
[0025]
步骤3-2中所述的将所述输入信号即离散时域信号样本直接降采样l倍的方法包括：
[0026]
xi(n)＝x
i-1
(l
·
n)
[0027]
其中，xi(n)表示迭代i次的离散时域信号样本，n表示所述离散时域信号样本的索引编号，下标i表示迭代索引，索引以i＝1开头，x0(n)为i＝1。
[0028]
步骤3-3中所述的更新窄带干扰频率的估计值的方法包括：
[0029]
在第i次迭代中，假设峰值k
nbi_hy
与上一次迭代中获得的峰值功率频点索引k
i-1,nbi
的关系如下：
[0030]knbi_hyp
＝l
·
κ
i-1,nbi
[0031]
将第i次迭代中获得的峰值功率频点索引k
i,nbi
更新为：
[0032][0033]
其中，w代表半搜索范围，
[0034]
g次迭代后，还原出真实窄带干扰频率k
nbi
的方法如下：
[0035]knbi
＝κ
i＝g,nbi
/lg。
[0036][0037]
步骤3-3中所述的迭代终止条件包括：
[0038]
监测κ
i,nbi
附近2个频点的离散傅立叶变换dft输出功率；当两个频点的功率都小于
阈值时，即：
[0039]
max(p(κ
i,nbi-1),p(κ
i,nbi
1))《p
α
[0040]
则表示达到迭代终止条件，即结束整个迭代估计的流程；其中，p(κ
i,nbi-1)表示离散傅立叶变换在κ
i,nbi-1这个频点的功率输出，p(κ
i,nbi
1)表示离散傅立叶变换在κ
i,nbi
1这个频点的功率输出，p(κ
i,nbi
)表示离散傅立叶变换dft在κ
i,nbi
这个频点的功率输出，p
α
为阈值。
[0041]
步骤3-3中所述的阈值p
α
预设为常数或p
α
＝β
·
p(κ
i,nbi
)；其中，β是预先设定的比率。
[0042]
本发明在上述步骤之外，还提出了一种替代步骤3中所述的根据所述频点索引计算得到窄带干扰频率的方法，包括：
[0043]
根据所述频点索引取得第一离散值即第一频点k0和第二离散值即第二频点k1，具体方法包括：
[0044][0045]
真实窄带干扰频率k
nbi
满足k0≤k
nbi
≤k1；
[0046]
计算第一频点k0的功率值p(k0)和第二频点k1的功率值p(k1)的比率∈：
[0047]
∈＝p(k0)/p(k1)，0《∈《∞
[0048]
根据功率值的比率∈估计窄带干扰频率k
nbi
，具体包括：
[0049]knbi
＝k0 φ
[0050]
其中，δ表示分数指数，也就是整数以外不满整数的数量，0≤δ《1。
[0051]
所述的δ的计算方法包括：
[0052]
根据反函数ω(∈)计算，或查表确定。
[0053]
所述的反函数ω(∈)，计算方法为：ω(∈)＝ψ-1
(δ)；即δ＝ω(∈)
[0054][0055]
有益效果：
[0056]
1、基于dft的方式且采用陷波滤波器处理nbi时，透过低成本的nbi频率估计方式，使陷波滤波器能在精确位置发挥作用，藉此提升通信效能，来保证在干扰下的通信可靠度。
[0057]
2、本发明提升了估算干扰频率位置的精准度，透过简约不繁杂的数学转换式，在dft上进行二分逼近法，并搭配合理的阀值与搜索范围，来有效地提升nbi频率估计精度。
附图说明
[0058]
下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。
[0059]
图1为本发明中窄带干扰频率估计的流程示意图。
[0060]
图2为窄带干扰与通信信号示意图。
[0061]
图3为相干信号的dft结果示意图。
[0062]
图4为非相干信号的dft以及包络示意图。
[0063]
图5为反函数示意图。
[0064]
图6为输入整数倍频点指标与输出混迭频点指标关系对照示意图。
[0065]
图7为每次迭代中更新的估计nbi频率示意图。
具体实施方式
[0066]
为了处理强nbi，一般来说使用陷波滤波方法来处理是个有效的手段。在实际实现中，可以考虑无限脉冲响应陷波滤波器(无限脉冲响应iir又称递归型滤波器)，它可在陷波频率附近提供深度衰减和尖锐的滤波器频率响应。
[0067]
由于nbi的强大功率，需要确保陷波滤波器的陷波频率与nbi频率足够接近，如此一来陷波滤波器才能够有效地消除/抑制nbi，同时最大限度地减少其他频率的有用信息的破坏与失真。也就是说，nbi的准确频率估计在评估过程中起着关键作用。这也是本发明的重点。
[0068]
一般而言，接收器应会安排一个时隙来侦测与识别nbi的存在。执行评估过程不仅检查是否存在接近单音(也就是nbi)的频率，还需估计出此一单音的频率。通常采用离散傅里叶变换(discrete fourier transform，缩写为dft)将时域波形转换为频域。并且，若有某个频率点(frequency bin)的功率远大于其他频率点，则代表检测到强干扰。通常，在这侦测nbi的过程中，计算所有频率点的平均功率，并将此一平均功率乘上一个预定的比率后作为检测阈值。如果发现存在着频率点的功率大于此一阈值，则识别为nbi。然后，则启动这个窄带干扰去除或者是抑制。总之，干扰信号的分析与确定对于nbi消除和避免技术都是非常需要的预处理过程。
[0069]
考虑dft的原因是dft通常已经在通信系统中存在着，特别是对于基于ofdm(orthogonal frequency division multiplexing，正交频分复用)的系统。此外，当分析的离散信号是周期性的并且与dft相干时，dft提供了相当准确的检测结果。然而，对于一个正常的信号，这个相干的条件通常不能满足。在这种情况下，dft频谱分析的结果并不精准，这种误差的主要影响被称为dft频谱泄漏。
[0070]
实施例：
[0071]
考虑幅度为a、频率f
nbi
和初始相位φ的单音频率。离散时域信号样本x(n)可以表示为
[0072]
x(n)＝a
·
cos(2π
·fnbi
·
(n
·
ts) φ)
[0073]
其中n为整数，ts为采样频率fs＝1/ts对应的采样周期。假设n是dft大小(即n点dft，n也是样本数)，如果满足如下条件，则不存在上述dft频谱泄漏问题：
[0074]
1)满足离散信号周期性条件；
[0075]
2)满足条件n
·
ts＝k
nbi
·
τ
nbi
；
[0076]
其中，nbi的频点频率k
nbi
必须是整数，指的是总周期t＝nts的离散时域信号样本x(n)的周期数，τ
nbi
＝1/f
nbi
是采样信号频率的周期；
[0077]
若满足条件，则输入信号与dft相干，而dft输出频谱仅由一个频率为f
nbi
(或简称为k
nbi
的整数)的频谱分量组成，如图3所示，为相干信号的dft结果示意图。在本例中，k
nbi
＝13。由于是相干信号，dft输出能量只出现在k
nbi
＝13这个频点上，其馀所有的频点能量均为零，也就是没有频谱泄漏。
[0078]
相反的，当不满足上述相干条件时，就会得到非相干信号。值得注意的是，k
nbi
将是采样率fs的“归一化”频率，因为k
nbi
＝n
·
(f
nbi
/fs)。在下文中，将f
nbi
和k
nbi
一起使用，因为它们都代表估计的nbi频率的频率。
[0079]
大多数情况下，在nbi频率的实际估计中存在不连贯的情况。也就是说，k
nbi
不是整数。信号的“真实”频率(频率为f
nbi
)位于dft频率轴的两个离散值(表示为第一离散值k0和第二离散值k1)之间，第一离散值k0和第二离散值k1表示如下：
[0080]
且
[0081]
输入信号的频率可以表示为k
nbi
＝k0 δ，其中频率位移δ在0《δ《1.的区间内。泄漏导致在频率k
nbi
周围等距分布的无限宽谐波分量频谱。它们的大小对应于包络函数
[0082]
sin(t)/t＝sinc(t),
[0083]
如图4所示，为非相干信号的dft以及包络示意图。不同于图3，非相干信号会在dft的输出中观测到许多带有能量的频点，也就是说，测试信号虽然是单频，却在许多dft输出频点中看到能量，这就是频谱泄漏的问题，无法由dft的输出精准地估测出单频信号的频率。因此，对于一个非相干信号使用dft获得的频谱不是原始信号的实际频谱，而是一个拖尾版本。看起来像是一个频率区间的信号能量泄漏到其他频率区间。在评估nbi频率时，频谱泄漏现象会导致在估算频率时出现问题。
[0084]
第一部分：
[0085]
在所有dft输出频点中找到输出最大功率的频率频点索引k
peak
，方法为：
[0086][0087]
其中，p(k)表示dft输出频点索引k的频点功率值，d表示频点索引k的查找范围；
[0088]
然后，频点索引即第一离散值k0和第二离散值k1为：
[0089][0090]
也就是说，第一频点k0的功率p(k0)和第二频点k1的功率p(k1)是所有dft输出频点中对应于2个最大功率的2个频点，并且k0《k1。事实上，nbi的频点频率为k
nbi
，它不一定是整数，因为nbi并不总是非相干信号。而且很容易得到这三个频率之间的关系为k0≤k
nbi
≤k1。
[0091]
计算获得k
nbi
的第一个建议解决方案是将p(k0)与p(k1)的功率进行比较。这2个功率值的比率∈表示为：
[0092]
∈＝p(k0)/p(k1)，0《∈《∞
[0093]
此一功率比率是分数指数δ的函数，且0≤δ《1和k
nbi
＝k0 δ。
[0094]
反函数ω(∈)为ω(∈)＝ψ-1
(δ)。以此为依据,使用功率比并且根据函数ω来计算δ的分数指数(因为δ＝ω(∈))。然后，就可以确定nbi频率为k
nbi
＝k0 δ。如图5所示，为此反函数ω(∈)输入∈和输出δ的关系。当∈为100＝1.0时，输出δ为0.5。当∈往无穷大靠近，输出的δ往1.0靠近。反之，当∈往0靠近时输出的δ则逐渐逼近0。在本发明中，进一步使用查找表来代替函数ω的计算。查找表的分辨率取决于nbi频率估计的精度要求。
[0095]
第二部分：
[0096]
采用一种迭代解决方案的方式，代替使用功率比∈作为估计nbi频率的工具。即使在噪声影响下执行nbi评估时，它也能提供更稳健的估计。
[0097]
具体来说，所述的迭代解决方案，如图1所示，包括以下步骤：
[0098]
步骤1，对于通信讯号，进行信号降采；
[0099]
步骤2，对降采后的信号，执行离散傅里叶转换；
[0100]
步骤4，窄带干扰频率更新，进入步骤5；
[0101]
步骤5，迭代终止检查，不满足终止条件的，进入步骤6；否则，迭代终止，完成所述基于离散傅立叶变换的低成本窄带干扰频率估计；
[0102]
步骤6，迭代数加1，进入步骤1。
[0103]
首先使用一个简单的例子来说明描述本发明的概念。另q[v]表示一个未知的实数值v的舍入(四舍五入)算子。例如q[v]＝3，我们只知道2.5≤v《3.5，但不知道v本身的值。两者之间的误差是|q[v]-v|≤0.5是舍入误差。基于同样的舍入算子q[v]，我们要够得到更好(更精确)估计v的一种方法是将v加倍后，再去进行舍入运算。也就是说，一个更好的估计表示为
[0104][0105]
这很直观，意味着估计误差减小了。
[0106]
在执行nbi频率估计的第一次迭代之前，在所有dft输出频点中找到承载最大功率的频率频点索引k
peak
。它表示为：
[0107]k0,nbi
＝k
peak
.
[0108]
下标“0”指的是搜索过程的迭代索引。对于每次迭代，将时域信号样本直接降采样2，也就是(在本实施例中采用l＝2)：
[0109]
xi(n)＝x
i-1
(2n),
[0110]
其中，xi(n)表示迭代i次的离散时域信号样本，下标i表示迭代索引，索引以i＝1开头，x0(n)为i＝1。众所周知，时域信号的降采样2倍，形同加倍nbi的频率。dft应用于具有相同dft点(样本)的时域信号xi(n)的抽取(降采样)版本。dft输出为p(k),k∈{0,1,
…
,n/2-1}。在迭代i中，最大功率频点将落在在假设峰值k
nbi_hyp
附近，k
nbi_hyp
与在先前迭代中获得的峰值功率频点索引κ
i-1,nbi
相关，
[0111]knbi_hyp
＝2
·
κ
i-1,nbi
[0112]
然后，将迭代i k
i,nbi
这个“整数倍nbi频率”更新为
[0113][0114]
这个倍增的频率是在对信号于时域中抽取采样的直接结果。随着迭代的进行，估计的频率指数κ
i,nbi
不断增长(由于在每次迭代中频率都持续地乘上2)，并且在数次迭代后κ
i,nbi
可能超过n/2。一旦超过n/2，就会出现频率混迭。为解决这个问题使用了而是dft输出p(k)的无限重复版本，它与p(k)的关系为
[0115][0116]
其中是频点索引k的混迭版本，|v|是v的绝对值，[v]指的是v的四舍五入。这个公式转换的目的在于将不断增大的频点索引k，折回到[0～n/2-1]的这个范围内。如此一来，我们也才能由dft的输出中找到对应的来进行处理。频率混迭映像方程也如图6所示。
[0117]
图6展示了频点索引k与其混迭版本之间的关系。
[0118]
当执行了g次迭代，相乘的nbi频率将被估计为κ
i＝g,nbi
。除以倍频因子2g，最终估计的nbi频率是
[0119]knbi
＝κ
i＝g,nbi
/2g.
[0120]
若考虑fs的采样率，最终估计的nbi频率为f
nbi
＝κ
nbi
·
(fs/n).。
[0121]
下面举例说明一个数字示例。若测试输入频率f
nbi
＝1.9654khz。考虑采样率fs＝1mhz和dft大小n＝256，归一化nbi频率为k
nbi
＝50.31424。如图7所示，显示了在每次迭代中更新的估计nbi频率。首先当频率被乘上2倍时,频率变成100.62848。根据规则，必须在k
0,nbi
＝50的2倍附近，也就是在p(100),p(101)与p(102)中去找出最大功率的频点，得到结果为k
1,nbi
＝101。接着进入下一次递回，频率再被乘上2倍,频率变成201.25696。根据规则，必须在κ
1,nbi
＝101的2倍附近，也就是在p(201),p(202)与p(203)中去找出最大功率的频点，得到结果为κ
1,nbi
＝201。依此规则持续更新nbi频率，直到整个递回处理结束为止。
[0122]
第三部分：
[0123]
随着估计迭代的进行，频率分辨率增加，因此，nbi频率估计也会逐渐的精准。迭代终止可以是事先预定好的或是动态决定的。所谓的预定方案是使用恒定的迭代计数，而动态方法动态检测是否在每次迭代中去决定在该次迭代后，判断是否就可以结束整个迭代估计的流程。
[0124]
当nbi频率估计得足够精确，相应的泄漏将被相应地抑制。因此，可以监测κ
i,nbi
附近2个频率区间的dft输出功率。当两个频点的功率都小于阈值时，即
[0125]
max(p(κ
i,nbi-1),p(κ
i,nbi
1))《p
α
[0126]
则结束整个迭代估计的流程,其中，p(κ
i,nbi-1)表示dft在κ
i,nbi-1这个频点的功率输出，p(κ
i,nbi
1)表示dft在κ
i,nbi
1这个频点的功率输出，阈值p
α
可以是常数或p
α
＝β
·
p(κ
i,nbi
)；其中，β是预先确定的比率，p(κ
i,nbi
)表示dft在κ
i,nbi
这个频点的功率输出。
[0127]
第四部分：
[0128]
对于每次迭代，可以通过将时域信号抽取l来将音调频率乘以l。本实施例中，l＝2。由于频率乘以一个更大的因子，因此可望减少迭代次数但是同时保有相似的估计结果精度。然而，这是以更大的搜索范围为代价的(对于κ
i,nbi
)。具体来说，κ
i,nbi
需要被修改为：
[0129][0130]
其中w代表半搜索范围，
[0131]
具体实现中，本技术提供计算机存储介质以及对应的数据处理单元，其中，该计算机存储介质能够存储计算机程序，所述计算机程序通过数据处理单元执行时可运行本发明提供的一种基于离散傅立叶变换的低成本窄带干扰频率估计方法的发明内容以及各实施例中的部分或全部步骤。所述的存储介质可为磁碟、光盘、只读存储记忆体(read-only memory，rom)或随机存储记忆体(random access memory，ram)等。
[0132]
本领域的技术人员可以清楚地了解到本发明实施例中的技术方案可借助计算机程序以及其对应的通用硬件平台的方式来实现。基于这样的理解，本发明实施例中的技术
方案本质上或者说对现有技术做出贡献的部分可以以计算机程序即软件产品的形式体现出来，该计算机程序软件产品可以存储在存储介质中，包括若干指令用以使得一台包含数据处理单元的设备(可以是个人计算机，服务器，单片机，muu或者网络设备等)执行本发明各个实施例或者实施例的某些部分所述的方法。
[0133]
本发明提供了一种基于离散傅立叶变换的低成本窄带干扰频率估计方法的思路及方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。