CN103634026B - 一种基于fpga的数字零中频自适应陷波方法 - Google Patents

Abstract

本发明公开一种基于FPGA新型数字零中频自适应陷波方法，包括以下步骤：步骤1：利用FFT变换将采样的中频信号变换到频域；步骤2：在频域内对中频信号进行频谱统计量分析和门限判决，判断NBI个数及其中心频率;步骤3：如果没有干扰，转到步骤7；如果有干扰，则通过正交数字下变频器对采样的中频信号进行正交解调，并利用低通滤波器只保留采样中频带通信号的正频率分量，变换获得基带信号；步骤4：对基带信号进行频谱搬移，将第一个干扰点搬到零频，通过高通滤波器滤除干扰频谱分量，从而抑制掉第一个干扰点；步骤5：重复步骤4的操作，直至全部抑制其他干扰点；步骤6：将干扰抑制处理后的信号变换为实信号；步骤7：由基带处理模块进行解调输出。

Description

一种基于FPGA的数字零中频自适应陷波方法

技术领域

本发明属于通信技术领域，涉及扩频通信系统中强窄带干扰抑制，具体涉及一种基于FPGA的数字零中频自适应陷波器。

背景技术

直扩通信系统(DDDS)中，发射典型的低截获概率信号，信号所占频谱宽，功率密度低，对多种类型的干扰具有较好的抑制性能，但实际中往往受到系统带宽的限制，使系统本身固有的扩频增益不足以抑制大功率的强窄带干扰NBI(Narrow-Band Interference)。故在解扩前需要借助信号处理技术在不增加信号带宽的条件下，对强窄带干扰进行抑制，提高系统的抗干扰能力。

目前扩频系统NBI抑制技术主要分为三种：时域干扰抑制技术，变换域干扰抑制技术，码辅助干扰抑制技术。关于这三种技术，郭黎利，殷复莲，卢满宏等人发表的文章《DSSS/CDMA系统窄带干扰抑制技术概述[J]》(电子学报.2009,10:2248-2257)中对其优缺点进行了详细的介绍，这三种技术的共同的特点是设计各种不同算法的陷波器(notch filter)，先对中频信号进行各种方式的变换处理，然后检测信号中的NBI，并使用陷波器消除信号中的特定频率NBI，最后对中频信号进行恢复补偿，以便于后续的信号基带处理和解调解扩。图1给出了扩频系统基于这三种技术的传统NBI抑制模型，由图可以看出，这种方法将干扰检测和干扰抑制紧密结合在一起，在NBI检测前需要对中频的原始信号x(t)做预处理，无论系统是否有NBI，对x(t)都有一定程度的破坏，虽然在后期采用各种方法来补偿恢复中频信号，但增加了系统设计的复杂度和硬件开销，工程上实现困难。

发明内容

因此，针对上述的问题，本发明提出一种基于FPGA新型数字零中频自适应陷波器，简化设计算法，降低系统设计的复杂度和硬件开销，并有效抑制扩频系统中的NBI，从而解决现有技术之不足。

为了解决上述技术问题，本发明所采用的设计思想是，将干扰检测和干扰抑制模块完全分离，检测模块的任务仅仅是采用传统的NBI检测技术判断扩频系统是否有NBI，如果没有NBI，直接输出原始中频信号x(t)；如果有NBI，复制x(t)信号，再送到干扰抑制模块对NBI陷波抑制后输出x'(t)。这样，由于不需要对干扰检测处理过的信号进行各种恢复补偿，因此很大程度上降低了系统实现的复杂性，提高了系统的可靠性。本发明采用上述新的扩频系统NBI抑制模型，结合可编程逻辑门阵列(FPGA)的数字零中频处理技术，设计了一种新型的FPGA数字零中频自适应陷波方法以及陷波器。该陷波方法首先对采样后的数字中频信号进行基于门限估计的频域NBI检测，并对NBI进行统计量分析，找到NBI的归一化中心频率；然后采用数字正交解调将中频信号搬移到零中频(基带)，并在基带将NBI通过频谱搬移到零频陷波滤除，再经过信号重构，恢复中频带通信号(实信号)，最后送到基带处理模块解调输出。同时给出了陷波器的两个主要模块：基于门限估计的频域NBI检测模块和数字零中频NBI抑制模块的具体实现方法，实现对直接序列扩频系统(DSSS)中NBI的自适应陷波抑制。

具体的，本发明的一种基于FPGA新型数字零中频自适应陷波方法，包括以下步骤：

步骤1：利用FFT变换将采样的中频信号变换到频域；

步骤2：在频域内对中频信号进行频谱统计量分析和门限判决，判断NBI个数及其中心频率；

步骤3：如果没有干扰，转到步骤7对采样的原始中频信号进行直接输出；如果有干扰，则进行干扰抑制，通过正交数字下变频器(DDC)对采样的中频信号进行正交解调，并利用低通滤波器只保留采样中频带通信号的正频率分量，变换获得基带信号；

步骤4：对基带信号进行频谱搬移，将第一个干扰点搬到零频，通过高通滤波器滤除干扰频谱分量，从而抑制掉第一个干扰点；

步骤5：如果干扰点多于一个，重复步骤4的操作，即将第一次或第二次干扰抑制处理后的信号再进行频谱搬移，将干扰点搬到零频，通过高通滤波器滤除干扰频谱分量，从而抑制掉其他干扰点(第二个或第三个干扰点等等)；

步骤6：将干扰抑制处理后的信号变换为实信号(中频带通信号)；

步骤7：由基带处理模块进行解调输出。

其中，所述步骤2中，门限判决是采用基于变换频域的NBI门限检测估计方法。现有技术中，对于干扰检测门限的设计，常用的方法有N-sigma算法、σ²最大似然估计法、中值滤波法等。由于N-sigma算法在工程上比较容易实现，因此，本发明采用该种算法，并对该算法的有关参量进行一定的工程近似计算。具体的，该门限判决包括如下步骤：

步骤21：建立信号模型，对于直扩通信系统，在发射端通过将信号与伪随机(Pseudo-Noise，PN)序列相乘，把信号扩展在一个宽的频谱上；在接收端，通过将相同的PN序列与接收信号相乘，从而恢复信号。接收的中频信号x(t)一般由3部分组成：

x(t)＝s(t)+i(t)+n(t)       (1)

其中s(t)是扩频信号，i(t)是窄带干扰，n(t)是加性高斯白噪声。

x(t)以PN码元速率采样后，可表示为：

x(k)＝s(k)+i(k)+n(k)         (2)

其中s(k)和n(k)在时间上不相关，而与i(k)具有相关性；

步骤22：NBI抑制首先需要对中频信号进行采样，对采样的数字信号通过1024点FFT对接收信号进行频谱分析，通过对频谱信息进行统计分析，找出干扰点，然后进行NBI抑制。其目标是从较宽频率范围的信号中去除窄带干扰。这些干扰占据了较少的频率点，并且其幅值高于FFT的噪声台阶。在没有窄带干扰的情况下扩频信号加上噪声是一个近似高斯分布。由于干扰叠加在信号上，使频谱分布发生了畸变。在没有受到干扰影响的部分，信号的频谱分布保持着高斯分布的形状；在受到影响的频率部分，频谱分布在期望信号频谱的上方。在没有干扰的情况下，大多数信号小于均值(μ)加两倍标准方差(σ)；当有一系列干扰或干扰的总能量增大时，平均值相对于高斯分布部分的均值上升了，并且标准方差也由于干扰的存在而增大。这种畸变要求选择一个更靠近均值的门限来滤除与干扰信号相关的频率部分。因此，本发明滤除门限的确定方法如下：首先将FFT数据按照每1024个点划分为不同的数据块，对每个1024点的FFT数据块，计算数据块中每个频率点幅度的分贝值，然后计算该数据块的标准方差σ和均值μ；根据标准方差σ从预先设定的加权集合选取适当的加权因子N，设定干扰检测门为:

Th＝μ+Nσ        (3)

式中，N为加权因子；

σ分5个水平σ₀，σ₁，σ₂，σ₃，σ₄，N也分为N₀，N₁，N₂，N₃，N₄，N的选择与σ有关。

$N=\left\{\begin{array}{c}{N}_0,\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_0\\ N_1,{\mathrm{\&sigma;}}_0<\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_1\\ N_2,{\mathrm{\&sigma;}}_1<\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_2\\ N_3,{\mathrm{\&sigma;}}_2<\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_3\\ N_4,{\mathrm{\&sigma;}}_3<\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_4\end{array}\right.---\left(4\right)$

预设的加权因子集可根据不同的信道环境进行调整，这样信号在经过FFT后根据均值和标准方差选择N来得到干扰检测门限。在N-sigma算法中，N的确定和选取准则是算法设计的关键。N值取得过大，则估计门限过高，干扰泄漏严重；反之，N值取得过小，则估计门限过低，会产生误判。

σ越大表明干扰越强，需要一个小的比例系数以保持门限在噪声的顶部并且与之接近。每个频率点逐个与门限进行比较，大于门限的区域被记录并查找其中心位置，中心位置就是归一化干扰的中心频率。由于五个候选比例系数(N₀～N₄)的取值是由FPGA编程决定的，这使得算法可以方便的应用于其它系统应用中。

步骤23：方差和均值的统计量工程近似计算：干扰检测分析为每个FFT数据块计算σ和μ，这些统计量用来计算门限。为了计算这些统计量，每个FFT数据块输出分贝值为10log₁₀(|X(k)|)，其中|X(k)|为采样后的中频信号X(k)的幅度，这里采用对数坐标，用来减小由扩频信号、热噪声和干扰混合而成的信号的频谱畸变，同时降低为接下来计算统计量所需要的|X(k)|的数值精度要求为了降低硬件实现的复杂度，对10log₁₀(|X(k)|)和统计量的计算进行了一些近似。经仿真和实际验证，由近似计算产生的精度下降不影响从扩频信号和热噪声混合信号中识别出干扰信号。采样后的中频信号X(k)的幅值|X(k)|近似计算如下：

|X(k)|≈[max(ReX(k),ImX(k))+min(ReX(k),ImX(k))/4]          (5)

其中Re是信号的实部，Im是信号的虚部，上述近似运算整体误差为0.6％，最大误差为11.6％(分别在π/4，3π/4，5π/4，7π/4处)。

对数量化的近似计算如下：

10log₁₀(|X(k)|)＝10log₁₀(2)·log₂(|X(k)|)≈3log₂(|X(k)|)       (6)

可见幅度计算仅通过移位和加法运算便可完成，而避免了对硬件资源敏感的乘法和开方运算，上述对数运算整体误差0.5％，最大误差2.9％。

每个FFT数据块的均值和标准方差的计算如下：

$\mathrm{\&mu;}=\underset{k=0}{\overset{N_P-1}{\mathrm{\&Sigma;}}}\frac{\left(10{\log }_{10}\left(\right|X\left(k\right)\left|\right)\right)}{N_P}---\left(7\right)$

${\mathrm{\&sigma;}}^2=\frac{1}{N_P}[\underset{k=0}{\overset{N_P-1}{\mathrm{\&Sigma;}}}{(10\&CenterDot;{\log }_{10}\left(\right|X\left(k\right)\left|\right))}^2-\frac{1}{N_P}\left(\underset{k=0}{\overset{N_P-1}{\mathrm{\&Sigma;}}}{(10\&CenterDot;{\log }_{10}\left(\right|X\left(k\right)\left|\right))}^2\right)---\left(8\right)$

其中，N_P为每个FFT数据块的频点数。根据上述等式，累加器累加N_P个频点的10log₁₀(|X(k)|)和(10log₁₀(|X(k)|))²，每个块处理完后，累加值用来计算均值μ和方差σ²，σ可由σ²开方得来。

进一步的，所述步骤3中的信号抑制包括以下步骤：

步骤31，构造基带信号：令中频带通信号的数学表达式为：

$\begin{array}{c}S\left(t\right)=a\left(t\right)\cos [2\mathrm{\&pi;}{f}_{0}t+\mathrm{\&phi;}\left(t\right)]\\ =\frac{1}{2}a\left(t\right)e^{j[2\mathrm{\&pi;}{f}_{0}t+\mathrm{\&phi;}\left(t\right)]}+\frac{1}{2}a\left(t\right)e^{-j[2\mathrm{\&pi;}{f}_{0}t+\mathrm{\&phi;}\left(t\right)]}\end{array}---\left(9\right)$

其中a(t)是信息序列，f₀是载波中心频率，为载波相位，那么频移滤波后的基带信号S_B(t)的数学表达式为：

式中S_BI(t)和S_BQ(t)分别是频移后中

频带通信号的同向分量和正交分量(I分量和Q分量)，是频移因子；

步骤32：对干扰点进行抑制：

将基带信号S_B(t)包含的干扰成分记为f_I,那么基带信号与频移因子e^-j2πfIt相乘后得到f_I被搬到零频的基带信号S_B ^*(t)，S_B ^*(t)经高通滤波后得到抑制掉干扰f_I的基带信号S_B'(t)；抑制掉干扰f_I的基带信号S_B'(t)的计算公式如下：

(S_I+jS_Q)×[cos(ω₀t)+jsin(ω₀t)]

＝[S_I×cos(ω₀t)-S_Q×sin(ω₀t)]+j[S_I×sin(ω₀t)+S_Q×cos(ω₀t)]       (11)

式中ω₀＝-2πf_I；

步骤33：中频带通信号重构：

干扰抑制后的信号是基带信号(复信号)，为了后续处理需将其转换为实信号，通过分析中频带通信号S(t)、解析信号S_A(t)、基带信号S_B(t)间的相互关系，可得到由基带信号重构实信号的数学表达式。解析信号S_A(t)可表达为：

$\begin{array}{c}{S}_A\left(t\right)=S\left(t\right)+j\hat{S}\left(t\right)\\ =S_B\left(t\right)e^{j2\mathrm{\&pi;}{f}_{c}t}\\ =[S_{\mathrm{BI}}\left(t\right)+j{S}_{\mathrm{BQ}}\left(t\right)][\cos \left(2\mathrm{\&pi;}{f}_{c}t\right)+j\sin \left(2\mathrm{\&pi;}{f}_{c}t\right)]\\ =[S_{\mathrm{BI}}\left(t\right)\cos \left(2\mathrm{\&pi;}{f}_{c}t\right)-S_{\mathrm{BQ}}\left(t\right)\sin \left(2\mathrm{\&pi;}{f}_{c}t\right)\\ +j[S_{\mathrm{BI}}\left(t\right)\sin \left(2\mathrm{\&pi;}{f}_{c}t\right)+S_{\mathrm{BQ}}\left(t\right)\cos \left(2\mathrm{\&pi;}{f}_{c}t\right)]\end{array}---\left(12\right)$

其中是中频带通信号S(t)的希尔伯特变换,f_c为中频带通信号的中心频率，中频带通信号的带宽为2B。比较(12)式等式两端，即得所需实信号为：

S(t)＝S_BI(t)cos(2πf_ct)-S_BQ(t)sin(2πf_ct)。

本发明的一种基于FPGA新型数字零中频自适应陷波器，包括模数转换单元(ADC)、FPGA单元、配置逻辑电路、存储器和用户接口，中频输入信号经过模数转换单元(ADC)发送至FPGA单元，存储器和用户接口均与配置逻辑电路双向通讯连接，配置逻辑电路的输出端接于FPGA单元的输入端；其中FPGA单元包括FFT(快速傅立叶变换)模块、基于频谱的干扰检测分析模块、干扰自适应陷波抑制模块和基带信号处理模块，FFT模块的输出端接于干扰检测分析模块的输入端，干扰检测分析模块的输出端接于干扰自适应陷波抑制模块的输入端和基带信号处理模块的输入端，干扰自适应陷波抑制模块的输出端还接于基带信号处理模块的输入端。

本发明提供了一种新颖的基于FPGA数字零中频输入信号自适应陷波方法，它借助变换域门限干扰检测技术和数字正交解调的思想，对中频输入信号进行实时的采样和频谱分析，并通过程序预置参数，选取适合的抑制门限，将信号搬移到零频进行滤波，NBI滤除后，再将基带信号恢复成中频带通信号进行解调解扩，实现对DDDS的NBI自适应抑制。本发明大大简化了设计算法，降低了系统设计的复杂度和硬件开销，并有效抑制了扩频系统中的NBI。

附图说明

图1给出了扩频系统基于现有技术的传统NBI抑制模型；

图2为本发明的扩频系统NBI抑制模型；

图3为本发明的FPGA数字零中频自适应陷波器结构；

图4为本发明的窄带干扰抑制流程；

图5为本发明的由中频带通信号构造基带信号的流程图；

图6a为本发明的中频带通信号构造基带信号频谱搬移过程(中频带通信号的频谱)；

图6b为本发明的中频带通信号构造基带信号频谱搬移过程(移频后的中频带通信号的频谱)；

图6c为本发明的中频带通信号构造基带信号频谱搬移过程(低通滤波器幅频特性)；

图6d为本发明的中频带通信号构造基带信号频谱搬移过程(移频后的中频带通信号经低通滤波生成基带信号)；

图7为本发明的干扰点抑制过程；

图8a为本发明的干扰抑制过程的频谱(信号存在干扰f_I)；

图8b为本发明的干扰抑制过程的频谱(将干扰f_I搬移到零频)；

图8c为本发明的干扰抑制过程的频谱(高通滤波器幅频特性)；

图8d为本发明的干扰抑制过程的频谱(经高通滤波后干扰被抑制)；

图9为本发明的由基带信号重构实信号的过程；

图10a为本发明中的基带信号S_B(t)的频谱；

图10b为本发明中的中频带通信号S(t)的频谱；

图10c为本发明中的解析信号S_A(t)的频谱；

图11为在没有中频信号输入情况下的背景噪声；

图12为抑制前的信号频谱的仿真曲线；

图13为经过陷波抑制后的信号频谱的仿真曲线。

具体实施方式

现结合附图和具体实施方式对本发明进一步说明。

本发明针对扩频系统中现有窄带干扰(NBI，Narrow-Band Interference)抑制方法的不足，提出了一种新型窄带干扰抑制的方法模型，并给出了这种模型下的数字零中频自适应陷波器的具体实现方法。其思想来源于通信系统中的正交解调，在结构和算法上不同于一般处理技术。仿真和实测验证表明，采用这种陷波(notch filter)方法的直接序列扩频系统(DSSS)能够对强窄带干扰具有较好的抑制作用。

本发明针对这种方法的不足，提出一种新型扩频系统抑制NBI的方法。其设计思想是将干扰检测和干扰抑制模块完全分离，检测模块的任务仅仅是采用传统的NBI检测技术判断扩频系统是否有NBI，如果没有NBI，直接输出中频输入信号x(t)；如果有NBI，复制x(t)信号，再送到干扰抑制模块对NBI陷波抑制后输出x'(t)。图2给出了这种新方法的实现模型，由于不需要对干扰检测处理过的信号进行各种恢复补偿，因此很大程度上降低了系统实现的复杂性，提高了系统的可靠性。

本发明采用这种新的扩频系统NBI抑制模型，结合可编程逻辑门阵列(FPGA)的数字零中频处理技术，设计了一种新型的FPGA数字零中频自适应陷波器，并给出了这种陷波器的两个主要模块：基于门限估计的频域NBI检测模块和数字零中频NBI抑制模块的具体实现方法，实现对直接序列扩频系统(DSSS)中NBI的自适应陷波抑制。

1系统模型

随着FPGA硬件和软件无线电技术的发展，扩频通信系统的中频数字化处理技术越来越成熟，本发明设计的数字零中频自适应陷波器，对采样后的数字中频输入信号进行基于门限估计的频域NBI检测，并对NBI进行统计量分析，找到NBI的归一化中心频率。采用数字正交解调将中频输入信号搬移到零中频(基带)，并在基带将NBI通过频谱搬移到零频陷波滤除，再经过信号重构，恢复中频带通信号(实信号)，最后送到基带处理模块解调输出。

1.1FPGA零中频数字自适应陷波器模型

基于FPGA的数字零中频自适应陷波器结构如图3所示。

数字零中频陷波器由模数转换器(ADC)、FPGA、配置逻辑电路、存储器和用户接口等电路组成。其中FPGA内部由快速傅立叶变换(FFT)模块、基于频谱的干扰检测分析模块、干扰自适应陷波抑制模块和基带信号处理模块等部分组成。实现各种功能的程序存放在存储器中，通过查表对FPGA进行动态配置来实对现NBI信号频谱的统计量分析，数字正交解调，NBI零中频陷波抑制，带通信号重构和基带解调输出等功能。

1.2干扰检测抑制方法

首先对采样的中频输入信号进行FFT变换到频域，在频域内对中频输入信号进行频谱统计量分析和门限判决，判断NBI个数及其中心频率。如果没有干扰，对采样的原始中频输入信号直接由基带处理模块进行解调输出。在有干扰的情况下，对采样的中频输入信号由正交数字下变频器(DDC)进行正交解调，并通过低通滤波器只保留采样中频带通信号的正频率分量，变换成基带信号(零中频输入信号)。

基带信号进行频谱搬移，将第一个干扰点搬到零频，通过高通滤波器滤除干扰频谱分量，从而抑制掉第一个干扰点。如果干扰点多于一个，将第一次或第二次干扰抑制处理后的信号再进行频谱搬移，将干扰点搬到零频，通过高通滤波器滤除干扰频谱分量，从而抑制掉第二个或第三个干扰点。最后将干扰抑制处理后的信号变换成实信号(中频带通信号)由基带处理模块进行解调输出。其实现流程如图4所示。

2基于门限判决的NBI检测模块设计

要实现对NBI的规避，对信道的准确估计是系统能够有效通信的基础。对电磁环境的检测可以采用变换域的功率谱估计来完成。谱估计技术经过多年的发展，己经有了许多成熟的方法，如周期图谱法，自回归(auto regressive，AR)法及小波/小波包谱估计技术等。通过将谱估计得到的功率谱与干扰检测门限进行比较，就可以确定干扰的位置。对于干扰检测门限的设计，常用的方法有N-sigma算法、σ²最大似然估计法、自适应多门限法、中值滤波法等。本发明设计采用工程上比较容易实现的N-sigma算法，并对该算法进行一定的工程应用改进。

2.1信号模型

本发明主要研究直扩通信系统(DDDS)，在发射端通过将信号与伪随机(Pseudo-Noise，PN)序列相乘，把信号扩展在一个宽的频谱上。在接收端，通过将相同的PN序列与接受信号相乘，恢复信号。DDDS接收的中频输入信号x(t)一般由3部分组成：

x(t)＝s(t)+i(t)+n(t)       (1)

其中s(t)是扩频信号；i(t)是窄带干扰；n(t)是加性高斯白噪声。

x(t)以PN chip速率采样后，可表示为：

x(k)＝s(k)+i(k)+n(k)        (2)

其中s(k)和n(k)在时间上不相关，而与i(k)具有相关性。

2.2检测方法

NBI抑制首先需要对中频输入信号进行采样，对采样的数字信号通过1024点FFT对接收信号进行频谱分析，通过对频谱信息进行统计分析，找出干扰点，然后进行NBI抑制。

NBI抑制的目标是从较宽频率范围的信号中去除窄带干扰。这些干扰占据了较少的频率点，并且其幅值高于FFT的噪声台阶。在没有窄带干扰的情况下扩频信号加上噪声是一个近似高斯分布。由于干扰叠加在信号上，使频谱分布发生了畸变。在没有受到干扰影响的部分，信号的频谱分布保持着高斯分布的形状；在受到影响的频率部分，频谱分布在期望信号频谱的上方。

在没有干扰的情况下，大多数信号小于均值μ加两倍标准方差σ；当有一系列干扰或干扰的总能量增大时，平均值相对于高斯分布部分的均值上升了，并且标准方差也由于干扰的存在而增大。这种畸变要求选择一个更靠近均值的门限来滤除与干扰信号相关的频率部分。

滤除门限的确定方法如下：首先对每个1024点的FFT数据块，计算块中每个频率点幅度的分贝值，然后计算该块的标准方差σ和均值μ。根据σ从预先设定的加权集合选取适当的加权因子N，设定干扰检测门为:

Th＝μ+Nσ(3)

式中:N为加权因子。

标准方差分5个水平σ₀，σ₁，σ₂，σ₃，σ₄，N也分为N₀，N₁，N₂，N₃，N₄，N的选择与σ有关。

$N=\left\{\begin{array}{c}{N}_0,\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_0\\ N_1,{\mathrm{\&sigma;}}_0<\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_1\\ N_2,{\mathrm{\&sigma;}}_1<\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_2\\ N_3,{\mathrm{\&sigma;}}_2<\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_3\\ N_4,{\mathrm{\&sigma;}}_3<\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_4\end{array}\right.---\left(4\right)$

预设的加权因子集可根据不同的信道环境进行调整，这样信号在经过FFT后根据均值和标准方差选择N来得到干扰检测门限。在N-sigma算法中，N的确定和选取准则是算法设计的关键。N值取得过大，则估计门限过高，干扰泄漏严重；反之，N值取得过小，则估计门限过低，会产生误判。

标准方差越大表明干扰越强，需要一个小的比例系数以保持门限在噪声的顶部并且与之接近。每个频率点逐个与门限进行比较，大于门限的区域被记录并查找其中心位置，中心位置就是归一化干扰的中心频率。由于五个候选比例系数(N₀～N₄)的取值是由FPGA编程决定的，这使得算法可以方便的应用于其它系统应用中。

2.3方差均值的统计量工程近似计算

干扰检测分析为每个FFT数据块计算均值和标准方差，这些统计量用来计算门限。为了计算这些统计量，每个FFT数据块输出分贝值为10log₁₀(|X(k)|)，其中|X(k)|为采样后的中频信号X(k)的幅度这里采用对数坐标，用来减小由扩频信号、热噪声和干扰混合而成的信号的频谱畸变，同时降低为接下来计算统计量所需要的|X(k)|的数值精度要求。为了降低硬件实现的复杂度，对10log₁₀(|X(k)|)和统计量的计算进行了一些近似。经仿真验证，由近似计算产生的精度下降不影响从扩频信号和热噪声混合信号中识别出干扰信号。幅度的近似计算如下：

|X(k)|≈[max(ReX(k),ImX(k))+min(ReX(k),ImX(k))/4]        (5)

其中Re是信号的实部，Im是信号的虚部，上述近似运算整体误差为0.6％，最大误差为11.6％(分别在π/4，3π/4，5π/4，7π/4处)。

第二个近似计算是对数量化，式子：

10log₁₀(|X(k)|)＝10log₁₀(2)·log₂(|X(k)|)≈3log₂(|X(k)|)        (6)

可见幅度计算仅通过移位和加法运算便可完成，而避免了对硬件资源敏感的乘法和开方运算，上述对数运算整体误差0.5％，最大误差2.9％。

每个FFT数据块的均值和标准方差的计算如下：

$\mathrm{\&mu;}=\underset{k=0}{\overset{N_P-1}{\mathrm{\&Sigma;}}}\frac{\left(10{\log }_{10}\left(\right|X\left(k\right)\left|\right)\right)}{N_P}---\left(7\right)$

${\mathrm{\&sigma;}}^2=\frac{1}{N_P}[\underset{k=0}{\overset{N_P-1}{\mathrm{\&Sigma;}}}{(10\&CenterDot;{\log }_{10}\left(\right|X\left(k\right)))}^2-\frac{1}{N_P}\left(\underset{k=0}{\overset{N_P-1}{\mathrm{\&Sigma;}}}{(10\&CenterDot;{\log }_{10}\left(\right|X\left(k\right)\left|\right))}^2\right)]---\left(8\right)$

其中N_P为每个FFT数据块的频点数，根据上述等式，累加器累加N_P个频点10log₁₀(|X(k)|)和(10log₁₀(|X(k)|))²，每个块处理完后，累加值用来计算均值μ和方差σ²，标准方差σ可由σ²开方得到。

3基于数字零中频的干扰抑制模块设计

为了滤除干扰点的信号，需要将采样的中频输入信号进行正交分解，分解之后中频带通信号被搬移到零中频，通过低通滤波得到基带信号，这些变换过程由FPGA中的数字正交下变频器DDC来完成。基带信号通过频谱搬移在零频陷波，滤除强窄带干扰信号，为了便于后续的解调解扩处理，需要将基带信号再重构成中频带通信号(实信号)。

3.1构造基带信号

本设计中系统接收中频输入信号为70MHz，数字采样速率为61.44MHz，经模数A/D数字化后中频的位置是F₀＝8.56MHz，故带通滤波器中心频率选择为8.56MHz，滤波器带宽则根据不同的信息速率进行选择，等于不扩频时窄带调制信号的带宽。为了滤除干扰点的信号，需要将中频带通信号变换成基带信号，通过DDC正交解调，采样信号与片内正余弦查找表(SINCOS LUT)输出的正交信号相乘(混频)，再通过低通滤波器得到基带信号，如图5所示。

设中频带通信号的数学表达式为：

$\begin{array}{c}S\left(t\right)=a\left(t\right)\cos [2\mathrm{\&pi;}{f}_{0}t+\mathrm{\&phi;}\left(t\right)]\\ =\frac{1}{2}a\left(t\right)e^{j[2\mathrm{\&pi;}{f}_{0}t+\mathrm{\&phi;}\left(t\right)]}+\frac{1}{2}a\left(t\right)e^{-j[2\mathrm{\&pi;}{f}_{0}t+\mathrm{\&phi;}\left(t\right)]}\end{array}---\left(9\right)$

其中a(t)是信息序列，f₀是载波中心频率，为载波相位，那么频移滤波

式中S_BI(t)和S_BQ(t)分别是频移后中频带通信号的同向分量和正交分量(I分量和Q分量)，是频移因子。。

由中频带通信号构造基带信号频谱搬移过程如图6所示，图6a为采样中频带通信号的频谱，图6b为移频后信号的频谱，图6c为低通滤波器幅频特性，图6d为移频信号经低通滤波生成基带信号的频谱。

其中中频带通信号的带宽为2B，f_C为其中心频率。

3.2干扰点抑制

假设基带信号S_B(t)包含干扰成分f_I，那么基带信号与频移因子相乘后得到f_I被搬到零频的基带信号S_B ^*(t)，S_B ^*(t)经高通滤波后得到抑制掉干扰f_I的基带信号S_B'(t)。该过程如图7所示。抑制掉干扰f_I的基带信号S_B'(t)的计算公式如下：

(S_BI+jS_BQ)×[cos(ω₀t)+jsin(ω₀t)]             (11)

＝[S_BI×cos(ω₀t)-S_BQ×sin(ω₀t)]+j[S_BI×sin(ω₀t)+S_BQ×cos(ω₀t)]

式中ω₀＝-2πf_I，干扰抑制过程中S_B(t)、S_B ^*(t)、S_B'(t)频谱如图8所示。

3.3中频带通信号重构

干扰抑制后的信号是基带信号(复信号)，为了后续处理需将其转换为中频带通信号(实信号)，通过分析中频带通信号S(t)、解析信号S_A(t)、基带信号S_B(t)间的相互关系，可得到由基带信号重构实信号的数学表达式。解析信号S_A(t)可表达为：

$\begin{array}{c}{S}_A\left(t\right)=S\left(t\right)+j\hat{S}\left(t\right)............\left(12\right)\\ =S_B\left(t\right)e^{j2\mathrm{\&pi;}{f}_{c}t}............\left(13\right)\\ =[S_{\mathrm{BI}}\left(t\right)+j{S}_{\mathrm{BQ}}\left(t\right)][\cos \left(2\mathrm{\&pi;}{f}_{c}t\right)+j\sin \left(2\mathrm{\&pi;}{f}_{c}t\right)]............\left(14\right)\\ =[S_{\mathrm{BI}}\left(t\right)\cos \left(2\mathrm{\&pi;}{f}_{c}t\right)-S_{\mathrm{BQ}}\left(t\right)\sin \left(2\mathrm{\&pi;}{f}_{c}t\right)\\ +j[S_{\mathrm{BI}}\left(t\right)\sin \left(2\mathrm{\&pi;}{f}_{c}t\right)+S_{\mathrm{BQ}}\left(t\right)\cos \left(2\mathrm{\&pi;}{f}_{c}t\right)]............\left(15\right)\end{array}$

令上式(12)与(15)实部相等，即得所需实信号为：

S(t)＝S_BI(t)cos(2πf_ct)-S_BQ(t)sin(2πf_ct)      (16)

其中是中频带通信号S(t)的希尔伯特变换,f_c为中频带通信号的中心频率。

由基带信号重构实信号的过程如图9所示。图10表示中频带通信号S(t)、解析信号S_A(t)、基带信号S_B(t)的频谱。

系统验证采用半实物仿真平台，模型采用数字微波接收机的数字中频部分。相关的仿真和测试参数如下：数据源采用均匀分布的随机信号，扩频码为平衡的1023位gold序列作，码速率为10.24M，信息调制方式是BPSK，扩频处理增益是20dB，接收中频信号为70MHz，数字样点采样速率为61.44MHz，经A/D数字化后中频的位置是8.56MHz，带通滤波器中心频率选择为8.56MHz，滤波器带宽等于不扩频时窄带调制信号的带宽。

仿真测试结果如下：在没有中频信号输入情况下的背景噪声见图11所示。在有中频信号输入的情况下，在陷波处理前，我们可以看到中频信号的频谱有两个明显的NBI尖峰，如图12所示。经过中频自适应陷波处理后，我们在信号频谱上可以看到NBI得到很好的抑制，如图13所示。

本发明重点研究了一种新颖的基于FPGA数字零中频输入信号自适应陷波技术，它借助变换域门限干扰检测技术和数字正交解调的思想，对中频输入信号进行实时的采样和频谱分析，并通过程序预置参数，选取适合的抑制门限，将信号搬移到零频进行滤波，NBI滤除后，再将基带信号恢复成中频带通信号进行解调解扩，实现对DDDS的NBI自适应抑制。文中给出了具体的自适应陷波的理论计算推导过程和相关频谱分析。通过仿真分析验证，证明该方法对NBI抑制的有效性。

尽管结合优选实施方案具体展示和介绍了本发明，但所属领域的技术人员应该明白，在不脱离所附权利要求书所限定的本发明的精神和范围内，在形式上和细节上可以对本发明做出各种变化，均为本发明的保护范围。

Claims (3)

1.一种基于FPGA新型数字零中频自适应陷波方法，包括以下步骤：

步骤1：利用FFT变换将采样的中频信号变换到频域；

步骤2：在频域内对中频信号进行频谱统计量分析和门限判决，判断强窄带干扰NBI个数及其中心频率；

步骤3：根据门限判决结果，如果没有干扰，转到步骤7对采样的原始中频信号进行直接输出；如果有干扰，则进行NBI干扰抑制，通过正交数字下变频器对采样的中频信号进行正交解调，并利用低通滤波器只保留采样中频带通信号的正频率分量，变换获得基带信号；

步骤4：对基带信号进行频谱搬移，将第一个干扰点搬到零频，通过高通滤波器滤除干扰频谱分量，从而抑制掉第一个干扰点；

步骤5：如果干扰点多于一个，重复步骤4的操作，即将第一次或第二次干扰抑制处理后的信号再进行频谱搬移，将干扰点搬到零频，通过高通滤波器滤除干扰频谱分量，从而抑制掉其他干扰点；

步骤6：将干扰抑制处理后的信号变换为实信号；

步骤7：由基带处理模块进行解调输出。

2.根据权利要求1所述的基于FPGA新型数字零中频自适应陷波方法，其特征在于：所述步骤2中的门限判决包括如下步骤：

步骤21：建立信号模型，对于直扩通信系统，在发射端通过将信号与伪随机序列相乘，把信号扩展在一个宽的频谱上；在接收端，通过将相同的PN序列与接收信号相乘，从而恢复信号；其中，接收的中频信号x(t)一般由3部分组成：

x(t)＝s(t)+i(t)+n(t)   (1)

其中s(t)是扩频信号，i(t)是窄带干扰，n(t)是加性高斯白噪声；

x(t)以PN码元速率采样后，表示为：

x(k)＝s(k)+i(k)+n(k)   (2)

其中s(k)和n(k)在时间上不相关，而s(k)与i(k)有相关性；

步骤22：将FFT数据按照每1024个点划分为不同的数据块，对每个1024点的FFT数据块，计算数据块中每个频率点幅度的分贝值，然后计算该数据块的标准方差σ和均值μ；根据标准方差σ从预先设定的加权集合选取适当的加权因子N，设定干扰检测门为:

Th＝μ+Nσ   (3)

式中，N为加权因子；

σ分5个水平σ₀，σ₁，σ₂，σ₃，σ₄，N分为N₀，N₁，N₂，N₃，N₄，N的选择与σ有关：

$N=\left\{\begin{array}{c}{N}_0,\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_0\\ N_1,{\mathrm{\&sigma;}}_0<\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_1\\ N_2,{\mathrm{\&sigma;}}_1<\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_2\\ N_3,{\mathrm{\&sigma;}}_2<\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_3\\ N_4,{\mathrm{\&sigma;}}_3<\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_4\end{array}\right.---\left(4\right)$

步骤23：对步骤22中的标准方差σ和均值μ的统计量进行工程近似计算：干扰检测分析为每个FFT数据块计算标准方差σ和均值μ，每个FFT数据块输出分贝值为10log₁₀(|X(k)|)；采样后的中频信号X(k)幅值|X(k)|的近似计算如下：

|X(k)|≈[max(ReX(k),ImX(k))+min(ReX(k),ImX(k))/4]   (5)

其中Re是信号的实部，Im是信号的虚部；

对数量化的近似计算如下：

10log₁₀(|X(k)|)＝10log₁₀(2)·log₂(|X(k)|)≈3log₂(|X(k)|)   (6)

每个FFT数据块的均值和标准方差的计算如下：

$\mathrm{\&mu;}=\underset{k=0}{\overset{N_P-1}{\mathrm{\&Sigma;}}}\frac{\left(10{\log }_{10}\left(\right|X\left(k\right)\left|\right)\right)}{N_P}---\left(7\right)$

${\mathrm{\&sigma;}}^2=\frac{1}{N_P}[\underset{k=0}{\overset{N_P-1}{\mathrm{\&Sigma;}}}{(10\&CenterDot;{\log }_{10}\left(\right|X\left(k\right)\left|\right))}^2-\frac{1}{N_P}\left(\underset{k=0}{\overset{N_P-1}{\mathrm{\&Sigma;}}}{(10\&CenterDot;{\log }_{10}\left(\right|X\left(k\right)\left|\right))}^2\right)]---\left(8\right)$

其中，N_P为每个FFT数据块的频点数；根据上述等式，累加器累加N_P个频点的10log₁₀(|X(k)|)和(10log₁₀(|X(k)|))²，每个块处理完后，累加值用来计算均值μ和方差σ²，标准方差σ由σ²开方得到。

3.根据权利要求1所述的基于FPGA新型数字零中频自适应陷波方法，其特征在于：所述步骤3中的NBI信号抑制包括以下步骤：

步骤31，构造基带信号：令中频带通信号的数学表达式为：

$\begin{array}{c}S\left(t\right)=a\left(t\right)\cos [2\mathrm{\&pi;}{f}_{0}t+\mathrm{\&phi;}\left(t\right)]\\ =\frac{1}{2}a\left(t\right)e^{j[2\mathrm{\&pi;}{f}_{0}t+\mathrm{\&phi;}\left(t\right)]}+\frac{1}{2}a\left(t\right)e^{-j[2\mathrm{\&pi;}{f}_{0}t+\mathrm{\&phi;}\left(t\right)]}\end{array}---\left(9\right)$

其中a(t)是信息序列，f₀是载波中心频率，为载波相位，那么频移滤波后的基带信号S_B(t)的数学表达式为：

式中S_BI(t)和S_BQ(t)是原中频信号的I分量和Q分量，是频移因子；

步骤32：对干扰点进行抑制：将基带信号S_B(t)包含的干扰成分记为f_I,那么基带信号与频移因子相乘后得到f_I被搬到零频的基带信号s*(t)，s*(t)经高通滤波后得到抑制掉干扰f_I的基带信号s'(t)；抑制掉干扰f_I的基带信号s'(t)的计算公式如下：

$\begin{array}{c}(S_{\mathrm{BI}}+j{S}_{\mathrm{BQ}})\&times;[\cos \left({\mathrm{\&omega;}}_{0}t\right)+j\sin \left({\mathrm{\&omega;}}_{0}t\right)]\\ =[S_{\mathrm{BI}}\&times;\cos \left({\mathrm{\&omega;}}_{0}t\right)-S_{\mathrm{BQ}}\&times;\sin \left({\mathrm{\&omega;}}_{0}t\right)]+j[S_{\mathrm{BI}}\&times;\sin \left({\mathrm{\&omega;}}_{0}t\right)+S_{\mathrm{BQ}}\&times;\cos \left({\mathrm{\&omega;}}_{0}t\right)]\end{array}---\left(11\right);$

式中ω₀＝-2πf_I；

步骤33：中频带通信号重构：干扰抑制后的信号是基带信号，为了后续处理需将其转换为实信号，通过分析中频带通信号S(t)、解析信号S_A(t)、基带信号S_B(t)间的相互关系，即得到由基带信号重构实信号的数学表达式；解析信号S_A(t)表达为：

$\begin{array}{c}{S}_A\left(t\right)=S\left(t\right)+j\hat{S}\left(t\right)\\ =S_B\left(t\right)e^{j2\mathrm{\&pi;}{f}_{c}t}\\ =[S_{\mathrm{BI}}\left(t\right)+j{S}_{\mathrm{BQ}}\left(t\right)][\cos \left(2\mathrm{\&pi;}{f}_{c}t\right)+j\sin \left(2\mathrm{\&pi;}{f}_{c}t\right)]\\ =[S_{\mathrm{BI}}\left(t\right)\cos \left(2\mathrm{\&pi;}{f}_{c}t\right)-S_{\mathrm{BQ}}\left(t\right)\sin \left(2\mathrm{\&pi;}{f}_{c}t\right)]\\ +j[S_{\mathrm{BI}}\left(t\right)\sin \left(2\mathrm{\&pi;}{f}_{c}t\right)+S_{\mathrm{BQ}}\left(t\right)\cos \left(2\mathrm{\&pi;}{f}_{c}t\right)]\end{array}---\left(12\right)$

其中是中频带通信号S(t)的希尔伯特变换,f_c为中频带通信号的中心频率；

比较(12)式等式两端，即得所需实信号为：

S(t)＝S_BI(t)cos(2πf_ct)-S_BQ(t)sin(2πf_ct)。