CN109639612A - 一种基于非线性最小二乘法的zpw-2000信号解调方法 - Google Patents

一种基于非线性最小二乘法的zpw-2000信号解调方法 Download PDF

Info

Publication number
CN109639612A
CN109639612A CN201811451854.6A CN201811451854A CN109639612A CN 109639612 A CN109639612 A CN 109639612A CN 201811451854 A CN201811451854 A CN 201811451854A CN 109639612 A CN109639612 A CN 109639612A
Authority
CN
China
Prior art keywords
frequency
value
carrier frequency
low
signal
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201811451854.6A
Other languages
English (en)
Other versions
CN109639612B (zh
Inventor
武晓春
戈壁舟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Lanzhou Jiaotong University
Original Assignee
Lanzhou Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Lanzhou Jiaotong University filed Critical Lanzhou Jiaotong University
Priority to CN201811451854.6A priority Critical patent/CN109639612B/zh
Publication of CN109639612A publication Critical patent/CN109639612A/zh
Application granted granted Critical
Publication of CN109639612B publication Critical patent/CN109639612B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04L—TRANSMISSION OF DIGITAL INFORMATION, e.g. TELEGRAPHIC COMMUNICATION
    • H04L27/00—Modulated-carrier systems
    • H04L27/10—Frequency-modulated carrier systems, i.e. using frequency-shift keying
    • H04L27/14—Demodulator circuits; Receiver circuits
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04L—TRANSMISSION OF DIGITAL INFORMATION, e.g. TELEGRAPHIC COMMUNICATION
    • H04L27/00—Modulated-carrier systems
    • H04L27/0014—Carrier regulation
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04L—TRANSMISSION OF DIGITAL INFORMATION, e.g. TELEGRAPHIC COMMUNICATION
    • H04L27/00—Modulated-carrier systems
    • H04L27/10—Frequency-modulated carrier systems, i.e. using frequency-shift keying

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Signal Processing (AREA)
  • Digital Transmission Methods That Use Modulated Carrier Waves (AREA)

Abstract

本发明公开了一种基于非线性最小二乘法的ZPW‑2000信号解调方法。该方法基于非线性最小二乘法法则，包含数据预处理、粗估计、网格搜索及精确搜索五个步骤。本发明具有较高的解调准确度。高信噪比环境下采样持续时间大于等于0.11s的ZPW‑2000移频信号均可被本发明正确解调。同时，本发明具备实时解调能力，较好的抗白噪声能力以及较好的抗单频干扰能力。

Description

一种基于非线性最小二乘法的ZPW-2000信号解调方法
技术领域
本发明属于铁路信号轨道电路领域，涉及一种基于非线性最小二乘法的 ZPW-2000信号解调方法，应用于普速及配备CTCS-2级列控系统的列车车载设备中ZPW-2000轨道电路移频信号的解调。
背景技术
ZPW-2000轨道电路作为一种地-车不间断通信设备，被广泛应用于我国既有线及客运专线，在ZPW-2000轨道电路区段，车载TCR(track circuit reader)天线不断感应以钢轨为传输媒介的ZPW-2000移频信号，车载设备利用感应电压采样值解调得到低频码序，从而获知前方区段占用情况，指导列车安全运行。
由于ZPW-2000移频信号(以下简称ZPW-2000信号)与列车牵引回流共用钢轨作为传输媒介，不平衡牵引电流干扰对ZPW-2000信号影响较为严重。同时，相关规范要求车载设备接收信息的应变时间不能超过2s。因此，ZPW-2000信号解调方法应兼备抗干扰能力和实时解调能力。
ZPW-2000信号为相位连续二进制频移键控(CPBFSK,continuous phasemodulated binary frequency shift keying)信号，其二进制符号”0”、”1”始终处于交替变化状态，使得其频谱形态基本恒定。考虑采样过程，ZPW-2000信号s(n)的傅立叶级数展开形式为：
式中，A0、f0、f1分别为s(n)的包络幅值、载频频率和低频频率；m＝f1/fΔ为调制指数；fΔ为频偏，恒为11Hz；Ts为采样间隔。式(1)描述了ZPW-2000信号理想频谱的形态，依据该理想频谱形态，发展出一系列基于信号频域的解调方法。
不平衡牵引电流干扰主要包括牵引谐波干扰和暂态冲击干扰。牵引谐波干扰谐波频率间隔为50Hz左右，因此在ZPW-2000信号有效频带[f0-40，f0+40](Hz) 内只可能存在一个谐波频率分量。因此本发明将牵引谐波干扰描述为单频干扰。暂态冲击干扰频率成分复杂，不失一般地，本发明将其假定为高斯白噪声。因此，在本发明实施例的仿真试验中，将抗白噪声能力和抗单频干扰能力作为解调方法抗干扰能力的评判指标。
现有机车信号设备通常采用基于FFT(fast Fourier transform)的周期图法(亦称频域分析法)来解调ZPW-2000信号，主要有：文献“移频轨道电路测试系统的设计”(魏学业,汪希时,丁正庭.铁道学报,18(05):67-72,1996)提出的结合欠采样、 FFT滤波及补零FFT的解调方法，文献“铁路移频信号处理方法研究”(胡幸江,黄文君,何伟挺,谭平.仪器仪表学报,33(08):1729-1734,2012)提出的结合欠采样和补零FFT的解调方法(以下简称“欠采样补零法”)，文献“铁路移频信号检测系统设计与实现”(杨帆,戴胜华,刘泽.电子测量与仪器学报,24(05):500-505,2010)提出的 Zoom-FFT结合重心法插值的解调方法(以下简称“Zoom-FFT插值法”)。由于 ZPW-2000信号各频率成分分布较集中，而周期图法存在频谱泄漏现象，能量较大的载频分量的泄漏值容易使附近能量较小的边频分量谱峰位置发生偏移，当采样持续时间较小时，这种偏移现象会更加严重，此时周期图法的解调准确度会大幅降低，解调正确性难以保障。当列车高速通过较短的道岔区段以及信号受到暂态冲击干扰影响时，有效信号的采样持续时间长度被缩短，周期图法的解调正确性会受到挑战。因此，通过采用新的解调方法，提升解调准确度，进而降低解调方法对采样持续时间长度的限制，对于改善ZPW-2000轨道电路这一地-车通信方式的可靠性具有重要意义。
发明内容
为解决现有技术存在的上述问题，本发明提出基于非线性最小二乘法的 ZPW-2000信号解调方法，其解调准确度优于目前常用的周期图法而且具备实时解调能力。
一种基于非线性最小二乘法的ZPW-2000信号解调方法，该方法的实现目的为寻求以下价值函数的最大值：
式中，分别为载频频率f0、低频频f1估计值；x＝[x(0)x(1)x(2)…x(N-1)]T为接收信号采样序列；N为采样序列长度，序列下标从0开始；[·]T为矩阵转置符号；Z为下文定义的ZPW-2000信号近似模型矩阵形式中的参数矩阵；
依据ZPW-2000信号的傅立叶级数展开形式，定义ZPW-2000信号s(n)的近似模型如下：
式中，L为近似模型中频率分量个数，本发明令L＝5；ωl，l＝1,2,...,L为近似模型中各频率分量的频率值，其定义为ωl＝2π(f0+(l-3)f1)Ts，Ts为接收信号采样间隔； Al，l＝1,2,...,L和l＝1,2,...,L分别为近似模型中各频率分量的幅值和相位；al， l＝1,2,...,L和bl，l＝1,2,...,L共同组成近似模型中各频率分量的系数，它们与Al及的关系为： ZPW-2000信号近似模型的矩阵形式定义如下：
式中，为近似模型序列；Z为参数矩阵，定义如下：
Z＝[c(ω1) c(ω2) … c(ωL) s(ω1) s(ω2) … s(ωL)] (5)
式中，c(ωl)，l＝1,2,...,L和s(ωl)，l＝1,2,…,L分别为近似模型中各频率分量的余弦、正弦序列，对于频率分量ωl，c(ωl)和s(ωl)的定义分别如下：
c(ωl)＝[1cos(ωl) cos(2ωl) … cos(ωl(N-1))]T (6)
s(ωl)＝[0sin(ωl) sin(2ωl) … sin(ωl(N-1))]T (7)
α为近似模型中各频率分量的系数向量，定义如下：
α＝[a1 a2 … aL -b1 -b2 … -bL]T (8)
为求解价值函数的最大值，基于非线性最小二乘法的ZPW-2000信号解调方法的步骤为：
A.计算解调过程中与接收信号不相关的参数，以全局变量的形式存入内存；
B.依次提取内存中针对粗估计的144种载频、低频配置方案对应的参数，采用高斯消元法求解信号在不同载频、低频配置方案下各频率分量幅值，并计算相应的价值函数值。依据各频率分量幅值特征，排除不符合ZPW-2000信号频谱特征的载频、低频配置方案，在剩余方案中，选取价值函数值最大者对应的载频、低频值作为粗估计结果；
C.提取内存中针对网格搜索的相关参数中与载频、低频粗估计结果相对应的参数；采用FFT快速构造网格搜索所需的中间变量，采用高斯消元法求解信号在不同载频、低频配置方案下各频率分量幅值，并计算相应的价值函数值。选取价值函数值最大者对应的载频、低频值作为网格搜索结果；
D.在网格搜索结果附近，采用二维二分搜索算法搜寻价值函数峰值位置，获得载频、低频精确估计值及对应的各频率分量幅值；
E.汇总估计结果，通过条件判断，判断接收信号是否为ZPW-2000信号以及信号是否受到干扰，输出载频、低频及提示信息。
本发明具有以下优点：本发明的解调准确度优于目前常用的周期图法；在高信噪比条件下，本发明可正确解调采样持续时间大于等于0.11s的ZPW-2000移频信号的低频、载频频率，降低了解调方法对采样持续时间的限制；本发明具备实时解调能力，当信号长度为8000时，本发明在MATLAB环境下的平均运行时间为0.32s；本发明的抗白噪声性能优于目前常用的周期图法；本发明给出了具有容错特性的ZPW-2000信号频谱特征检测方法，当信号频谱形态被干扰信号破坏时，可发现译码错误，当信号受信干比大于等于1:3.6的单频干扰影响时，可正确解调。
本发明基于式(2)所示目标函数提出，详细给出了一种实现方法，使得解调方法能够满足实时解调需求。凡是以式(2)为目标函数的ZPW-2000移频信号解调方法都属于本发明的范围。
附图说明
图1为本发明流程图；
图2为本发明实施例1的仿真结果，即在采样持续时间为0.08s，以及以0.01s 为间隔从0.10s枚举至0.19s，以及以0.1s为间隔从0.2s枚举至1s共20种情形下，采用本发明方法，欠采样补零法及Zoom-FFT插值法解调ZPW-2000仿真信号所得的载频估计平均绝对误差结果；其中，每种采样持续时间下共有4320个高斯白噪声信噪比为30dB，采样频率为8000Hz的ZPW-2000仿真信号，涵盖了144种载频、低频配置方案，且对每种方案进行了30次MonteCarlo模拟；因此各解调方法的平均绝对误差为4320个绝对误差的平均值；
图3为本发明实施例1的仿真结果，与图2相对应，是实施例1仿真过程得出的低频估计平均绝对误差；
图4为本发明实施例1的仿真结果，与图2相对应，是实施例1仿真过程得出的错误译码次数统计数据；
图5为本发明实施例2的仿真结果，即当采样持续时间以0.1s为间隔从0.1s 枚举至1s时，采用本发明方法，欠采样补零法及Zoom-FFT插值法解调1000次载频、低频、采样频率分别为2001.4Hz，29Hz，8000Hz的ZPW-2000移频信号所花的平均运行时间。因此任意采样持续时间下各解调方法的平均运行时间为 1000个运行时间的平均值；
图6为本发明实施例3的仿真结果，即当高斯白噪声信噪比以2dB为间隔从-20dB枚举至50dB时，采用本发明方法，欠采样补零法及Zoom-FFT插值法解调ZPW-2000仿真信号所得的载频估计平均绝对误差结果。其中，每种高斯白噪声信噪比下共有4320个采样持续时间为0.3s，采样频率为8000Hz的ZPW-2000 仿真信号，涵盖了144种载频、低频配置方案，且对每种方案进行了30次Monte Carlo模拟。因此任意高斯白噪声信噪比下各解调方法的平均绝对误差为4320个绝对误差的平均值；
图7为本发明实施例3的仿真结果，与图6相对应，是实施例3仿真过程得出的低频估计平均绝对误差；
图8为本发明实施例3的仿真结果，与图6相对应，是实施例3仿真过程得出的错误译码次数统计数据；
图9为本发明实施例4的仿真结果，即当单频干扰的频率fsfi以1Hz为间隔从(2001.4-40)Hz枚举至(2001.4+40)Hz，幅值Asfi以0.2为间隔从0.2枚举至7时，采用本发明方法，欠采样补零法及Zoom-FFT插值法解调ZPW-2000信号与单频干扰叠加而成的仿真信号所得的译码错误次数统计结果。其中，ZPW-2000信号的采样频率、载频、低频、采样持续时间、高斯白噪声信噪比分别为8000Hz、 2001.4Hz、29Hz、1s、5dB，且在每种单频干扰情况下进行了30次Monte Carlo 模拟。因此任意频率、幅值的单频干扰下各解调方法的译码错误次数为30次仿真结果的统计值。
具体实施方式
下面对本发明的技术方案进行详细描述。
参考图1，一种基于非线性最小二乘法的ZPW-2000信号解调方法，该方法的实现目的为寻求以下价值函数的最大值：
式中，分别为载频频率f0、低频频f1估计值；x＝[x(0) x(1) x(2) … x(N-1)]T为接收信号采样序列；N为采样序列长度，序列下标从0开始；[·]T为矩阵转置符号；Z为下文定义的ZPW-2000信号近似模型矩阵形式中的参数矩阵。
由于ZPW-2000信号能量主要集中于载频及一、二次边频分量，依据ZPW-2000信号的傅立叶级数展开形式，定义如下ZPW-2000信号s(n)的近似模型如下：
式中，L为近似模型中频率分量个数，本发明令L＝5；ωl，l＝1,2,...,L为近似模型中各频率分量的频率值，其定义为ωl＝2π(f0+(l-3)f1)Ts，Ts为接收信号采样间隔； Al，l＝1,2,...,L和l＝1,2,...,L分别为近似模型中各频率分量的幅值和相位；al， l＝1,2,...,L和bl，l＝1,2,...,L共同组成近似模型中各频率分量的系数，它们与Al及的关系为： ZPW-2000信号近似模型的矩阵形式定义如下：
式中，为近似模型序列；Z为参数矩阵，定义如下：
Z＝[c(ω1) c(ω2) … c(ωL) s(ω1) s(ω2) … s(ωL)] (5)
式中，c(ωl)，l＝1,2,...,L和s(ωl)，l＝1,2,…,L分别为近似模型中各频率分量的余弦、正弦序列，对于频率分量ωl，c(ωl)和s(ωl)的定义分别如下：
c(ωl)＝[1cos(ωl) cos(2ωl) … cos(ωl(N-1))]T (6)
s(ωl)＝[0sin(ωl) sin(2ωl) … sin(ωl(N-1))]T (7)
α为近似模型中各频率分量的系数向量，定义如下：
α＝[a1 a2 … aL -b1 -b2 … -bL]T (8)
为求解价值函数的最大值，基于非线性最小二乘法的ZPW-2000信号解调方法的步骤为：
A.计算解调过程中与接收信号不相关的参数，以全局变量的形式存入内存。
B.依次提取内存中针对粗估计的144种载频、低频配置方案对应的参数，采用高斯消元法求解信号在不同载频、低频配置方案下各频率分量幅值，并计算相应的价值函数值。依据各频率分量幅值特征，排除不符合ZPW-2000信号频谱特征的载频、低频配置方案，在剩余方案中，选取价值函数值最大者对应的载频、低频值作为粗估计结果。
C.提取内存中针对网格搜索的相关参数中与载频、低频粗估计结果相对应的参数。采用FFT快速构造网格搜索所需的中间变量，采用高斯消元法求解信号在不同载频、低频配置方案下各频率分量幅值，并计算相应的价值函数值。选取价值函数值最大者对应的载频、低频值作为网格搜索结果。
D.在网格搜索结果附近，采用二维二分搜索算法搜寻价值函数峰值位置，获得载频、低频精确估计值及对应的各频率分量幅值。
E.汇总估计结果，通过条件判断，判断接收信号是否为ZPW-2000信号以及信号是否受到干扰，输出载频、低频及提示信息。
以下是基于非线性最小二乘法的ZPW-2000信号解调方法的详细实施步骤：
A.计算与接收信号不相关参数存入内存
(A.1)首先确认接收信号序列的长度N和采样频率fs。定义载频序列以及低频序列依次枚举载频和低频，形成144组编号为(p,q)的载频、低频配置方案分别构造对应的矩阵，·(p,q)表示编号为(p,q)的载频、低频配置方案下的数据。矩阵的参数配置方法参考式(3)和式(5)，其中的第一行元素采用CORDIC算法生成，的剩余元素由三角函数和角公式求得。由可得将144组与以全局变量的形式存入内存，供后续步骤(B.1)使用。
(A.2)定义载频序列2601.4,2598.7)(Hz)以及低频序列依次枚举载频和低频，形成144组编号为(p,q)的载频、低频配置方案对于编号(p,q) 的配置方案，载频和低频的网格搜索间隔均为ΔωAII，定义如下：
式中：
式中，为向下取整符号，意为四舍五入取整；以单位为Hz的实际频率值来描述，网格搜索中载频搜索频率值定义在(Hz)区间内，低频搜索频率值定义在(Hz)区间内；转化为数字频率，载频搜索频率值的定义如下：
式中，为：
当载频的网格搜索频率值确定后，即可以j为系数设定低频搜索值。j以及低频搜索值确定后各频率分量频率值l＝1,2,…,L，L＝5的定义如下所示：
式中，·(p,q,i,j)表示编号(p,q)的载频、低频配置方案下编号(i,j)的网格搜索值对应的数据；和的定义如下：
将式(16)至式(20)代入式(5)，可得矩阵计算出因此，144组载频、低频配置方案中，每组方案对应一组 M(p,q)，同时还包括组按次序以全局变量的形式将这些参数存入内存，供后续步骤(C.1)使用。
B.对ZPW-2000信号的载频、低频进行粗估计
(B.1)令JBImax＝0，f0c＝-3000(Hz)，f1c＝-3000(Hz)。确认长度为N的接收信号序列x＝[x(0)x(1)…x(N-1)]T。依次提取内存中(A.1)步骤存入的144组载频、低频配置方案对应的数据。其中，对于载频、低频配置方案其包含和采用高斯消元法求解方程：
将解得的αBI代入下式：
式中，Al，l＝2,3,4为各频率分量幅值，参考式(3)和式(8)由αBI求得；Athreshold为ZPW-2000信号包络幅值最小值，可实测获得；Ae为高信噪比下各频率分量幅值估计的最大误差，可实测获得，依据仿真结果，令Ae＝0.1；m为调制指数，其定义为r01为载频分量幅值与一次边频分量幅值比值的绝对值，定义如下：
若式(24)结果为真，则将αBI代入式(26)，得到该载频、低频配置方案的价值函数值：
否则令JBI＝0；此时，若JBI>JBImax，则更新粗估计价值函数及载频、低频粗估计值，即当完成144组载频、低频配置方案的枚举后，若f0c＝-3000(Hz)或f1c＝-3000(Hz)，则令并转至步骤E，否则转至步骤C。
C.采用网格搜索获得更精确的低频、载频估计值
(C.1)提取内存中(A.2)存入的数据中满足条件的载频、低频配置方案对应的计算：
ωshift＝ω0c-k0ΔωCI (27)
式中k0及其定义式中ω0c和ΔωCI定义分别如下：
定义向量xshift：
将序列x·xshift补零至M点，执行M点FFT得XM＝[XM(0)XM(1)XM(2)… XM(M-1)]。其中，·为向量点积符号。令JCImax＝0。依次提取(A.2)存入数据中载频、低频配置方案组以(i,j)为编号的网格搜索数据。对于其中某一网格搜索编号(i,j)，计算η＝[η1 η2 η3 ... η2L]T，其中：
式中，Re(·),Im(·)分别为取实部符号和取虚部符号；ξl，l＝1,2,…,L，L＝5的定义如下：
ξ3＝XM(k0+i) (35)
采用高斯消元法计算方程：
将αCI代入下式：
JCI＝ηTαCI (39)
此时，若JCI>JCImax，则JCImax＝JCI，ω0g＝ω0c+iΔωCI，ω1g＝(k1 (p,q)+j)ΔωCI，Δωg＝ΔωCI。当完成(2Lh+1)2组网格搜索值的枚举，转至步骤D。
D.采用精确搜索获得ZPW-2000信号载频、低频精确估计值
(D.1)精确搜索由二维二分搜索算法实现。定义载频频率二分搜索初始区间的左边界ω0l、右边界ω0r及区间半长Δω0定义为：ω0l＝ω0g-Δωg，ω0r＝ω0g+Δωg，Δω0＝Δωg。搜索终止条件为2Δω0≤ωrp，ωrp为解调方法的精度要求。载频频率二分搜索步骤如下：
Step 1：以载频频率ω0l、ω0r为参数进行低频频率二分搜索，分别得到对应的价值函数值J0l、J0r。令ω0m＝(ω0l+ω0r)/2。
Step 2：以载频频率ω0m为参数进行低频频率二分搜索，得到对应的低频估计值ω0m_1以及对应的频率分量系数α0m和价值函数值J0m。令Δω0＝Δω0/2，ω0m_last＝ω0m。若J0l<J0r，则J0l＝J0m，ω0m＝ω0m+Δω0，否则J0r＝J0m，ω0m＝ω0m-Δω0。
Step 3：若2Δω0≤ωrp，则结束二分搜索，进而载频、低频及各频率分量系数的精确估计值分别为否则回到 step 2。
低频频率二分搜索步骤与载频频率二分搜索基本一致。相对应地，低频频率二分搜索初始区间的左边界ω1l、右边界ω1r及区间半长Δω1定义为：ω1l＝ω1g-Δωg，ω1r＝ω1g+Δωg，Δω1＝Δωg。搜索终止条件为2Δω1≤ωrp。令代入的载频频率参数值为ω0p。低频频率二分搜索步骤如下：
Step 1：以ω0p为载频频率参数，ω1l、ω1r为低频频率参数，代入式(3)和式(5)，得到矩阵Z1l、Z1r，参考式(23)和式(26)，结合高斯消元法，得到对应的价值函数值J1l、J1r。令ω1m＝(ω1l+ω1r)/2。
Step 2：同样地，参考式(3)、式(5)、式(23)和式(26)求得载频频率ω0p，低频频率ω1m对应的价值函数值J1m以及频率分量系数α1m。令Δω1＝Δω1/2，ω1m_last＝ω1m。若J1l<J1r，则J1l＝J1m，ω1m＝ω1m+Δω1，否则J1r＝J1m，ω1m＝ω1m-Δω1。
Step 3：若2Δω1≤ωrp，则结束低频频率二分搜索，低频频率二分搜索输出的低频估计值，价值函数值以及频率分量系数分别为ω1m_last，J1m，α1m，否则回到 step 2。
E.汇总估计结果，通过条件判断，决定输出信息
(E.1)若或则提示信号接收错误，结束解调过程，否则将代入下式：
式中，l＝2,3,4为各频率分量幅值，参考式(3)和式(8)由求得；为载频分量幅值与一次边频分量幅值的比值，参考式(24)和式(25)求得。若式(40)结果为真，则正常显示载频、低频信息结束解调过程，否则显示载频、低频信息的同时，提示信号受到干扰，结束解调过程。
以下从具体仿真实验进一步说明本发明及其效果。
实施例1，为测试本发明的解调准确度，本实施例中，仿真生成了不同采样持续时间，144种载频、低频配置方案下，高斯白噪声信噪比、采样频率分别为 30dB、8000Hz的ZPW-2000仿真信号。其中，高斯白噪声信噪比的定义为σ2为高斯白噪声的方差。采用本发明方法，欠采样补零法及 Zoom-FFT插值法解调仿真信号，得到载频、低频平均绝对误差分别如图2、图 3所示。需指出，实施例1以及以下各实施例中，本发明的ωrp＝2π·(0.0005/8000)(rad/sample)，同时，通过补零，欠采样补零法及Zoom-FFT插值法中FFT的频域采样间隔为0.0005Hz。可看出，本发明方法的载频、低频估计平均绝对误差均小于欠采样补零法及Zoom-FFT插值法。当采样持续时间为0.3s 时，本发明的载频、低频平均绝对误差分别可低至0.0039Hz及0.0035Hz。因此，本发明的解调准确度优于目前常用的周期图法。
上述仿真过程中，若解调方法输出的载频或低频频率估计值绝对误差超过0.4Hz，则判定为译码错误。通过对每次仿真结果进行统计，在不同采样持续时间下，本发明方法，欠采样补零法及Zoom-FFT插值法译码错误次数统计结果如图4所示。当采样持续时间大于等于0.11s时，本发明译码错误次数可降为0次。相比之下，欠采样补零法及Zoom-FFT插值法都需要0.19s的采样持续时间才能保证信号正确解调。因此，本发明方法显著降低了解调方法对采样持续时间的限制。
实施例2，为测试本发明的实时性能，本实施例中，仿真生成了不同采样持续时间下，载频、低频、采样频率分别为2001.4Hz，29Hz，8000Hz的ZPW-2000 移频信号。采用本发明方法，欠采样补零法及Zoom-FFT插值法分别对每种仿真信号解调1000次，得到平均运行时间如图5所示。需指出，由于步骤A可在接收信号前完成，仿真中记录的运行时间不包括步骤A所花时间。可看出，虽然本发明的平均运行时间高于欠采样补零法及Zoom-FFT插值法，但未超过0.32s。因此本发明能使解调过程的应变时间不大于2s，基本可满足实时解调需求。由于本发明的粗估计、网格搜索过程中计算过程相互独立，利用相关硬件的并行计算功能，可进一步减少运行时间。此外，为了在满足实时解调需求的同时在时间上避开暂态冲击干扰的影响，可首先利用粗估计过程搜寻出粗估计结果能够保持一致的不重叠区间，之后再采用网格搜索及精确搜索获得精确的载频、低频估计值，进而可提高解调方法的可靠性。
实施例3，为测试本发明的抗白噪声能力，本实施例中，仿真生成了不同高斯白噪声信噪比，144种载频、低频配置方案下，采样持续时间、采样频率分别为0.3s、8000Hz的ZPW-2000移频信号。采用本发明，欠采样补零法及Zoom-FFT 插值法解调仿真信号，得到载频、低频平均绝对误差分别如图6、图7所示。可看出，相比于Zoom-FFT插值法，本发明需在信噪比大于等于-2dB后，载频、低频估计平均绝对误差才趋于稳定。这是因为本发明中，式(25)的判决条件可以甄别出部分译码错误情况，输出-3000Hz的载频、低频估计结果，因此部分译码错误情况不会被求平均过程掩盖，因此在信噪比小于-2dB时，本发明的载频、低频平均绝对误差高于Zoom-FFT插值法。
上述仿真过程中，若解调方法输出的载频或低频频率估计值绝对误差超过0.4Hz，则判定为译码错误。通过对每次仿真结果进行统计，在不同高斯白噪声信噪比下，本发明，欠采样补零法及Zoom-FFT插值法译码错误次数统计结果如图8所示。可看出，本发明在低信噪比下错误译码次数要小于前两种方法，在信噪比大于等于2dB时，本发明错误译码次数降为0。相比之下，欠采样补零法及 Zoom-FFT插值法分别需在信噪比大于等于16dB、4dB时方可保证正确解调。因此，本发明的抗白噪声能力优于欠采样补零法及Zoom-FFT插值法。同时，由图 6、图7与图8相比较可知，本发明中式(25)的判决条件尚不能甄别出所有的译码错误情况，因此当白噪声干扰较严重时，本发明的解调结果不再可信。
实施例4，由于29Hz的ZPW-2000信号边频幅值最小，最易受单频干扰影响，为测试本发明的抗单频干扰能力，本实施例仿真生成了低频频率为29Hz的 ZPW-2000信号以及幅值、频率分别为Asfi、fsfi的单频正弦干扰。两者相叠加，采用本发明，欠采样补零法及Zoom-FFT插值法解调该信号。当解调方法输出的载频或低频频率估计值绝对误差超过0.4Hz，则判定为译码错误。由此得到在不同幅值、频率单频干扰下各方法的译码错误次数统计情况如图9所示。图中，每个像素点颜色越深代表译码错误次数越多。可看出，当Asfi小于等于3.6时，本发明的译码错误次数均为0，可保证信号正确解调。因此，当信号受信干比大于等于1:3.6的单频干扰影响时，本发明可正确解调信号。
本发明是以文献"Fast fundamental frequency estimation:Making astatistically efficient estimator computationally efficient"(Nielsen J K,Jensen T L,Jensen J R,et al.Signal Processing,135:188-197,2017)所提的非线性最小二乘法为基础，结合 ZPW-2000移频信号特点提出一系列改进措施，而最终完成。此外，本发明参考了文献“A survey of CORDIC algorithms for FPGA based computers”(Andraka R. Proceedings of the 1998ACM/SIGDA sixth international symposium onField programmable gate arrays:191-200,1998)提出的CORDIC算法，用于计算任意角度的sin函数和cos函数值，以及文献“Linear algebra and its application 4th edition”(G. Strang.Cengage Learning:1-64,2006)介绍的高斯消元法，用于求解线性方程组。

Claims (7)

1.一种基于非线性最小二乘法的ZPW-2000信号解调方法，其特征在于，该方法的实现目的为寻求以下价值函数的最大值：
式中，分别为载频频率f0、低频频f1估计值；x＝[x(0) x(1) x(2) … x(N-1)]T为接收信号采样序列；N为采样序列长度，序列下标从0开始；[·]T为矩阵转置符号；Z为下文定义的ZPW-2000信号近似模型矩阵形式中的参数矩阵；
依据ZPW-2000信号的傅立叶级数展开形式，定义ZPW-2000信号s(n)的近似模型如下：
式中，L为近似模型中频率分量个数，本发明令L＝5；ωl，l＝1,2,...,L为近似模型中各频率分量的频率值，其定义为ωl＝2π(f0+(l-3)f1)Ts，Ts为接收信号采样间隔；Al，l＝1,2,...,L和分别为近似模型中各频率分量的幅值和相位；al，l＝1,2,...,L和bl，l＝1,2,...,L共同组成近似模型中各频率分量的系数，它们与Al及的关系为：ZPW-2000信号近似模型的矩阵形式定义如下：
式中，为近似模型序列；Z为参数矩阵，定义如下：
Z＝[c(ω1) c(ω2) … c(ωL) s(ω1) s(ω2) … s(ωL)] (5)
式中，c(ωl)，l＝1,2,...,L和s(ωl)，l＝1,2,…,L分别为近似模型中各频率分量的余弦、正弦序列，对于频率分量ωl，c(ωl)和s(ωl)的定义分别如下：
c(ωl)＝[1 cos(ωl) cos(2ωl) … cos(ωl(N-1))]T (6)
s(ωl)＝[0 sin(ωl) sin(2ωl) … sin(ωl(N-1))]T (7)
α为近似模型中各频率分量的系数向量，定义如下：
α＝[a1 a2 … aL -b1 -b2 … -bL]T (8) 。
2.根据权利要求1所述的一种基于非线性最小二乘法的ZPW-2000信号解调方法，其特征在于，为求解价值函数的最大值，基于非线性最小二乘法的ZPW-2000信号解调方法的步骤为：
A.计算解调过程中与接收信号不相关的参数，以全局变量的形式存入内存；
B.依次提取内存中针对粗估计的144种载频、低频配置方案对应的参数，采用高斯消元法求解信号在不同载频、低频配置方案下各频率分量幅值，并计算相应的价值函数值。依据各频率分量幅值特征，排除不符合ZPW-2000信号频谱特征的载频、低频配置方案，在剩余方案中，选取价值函数值最大者对应的载频、低频值作为粗估计结果；
C.提取内存中针对网格搜索的相关参数中与载频、低频粗估计结果相对应的参数；采用FFT快速构造网格搜索所需的中间变量，采用高斯消元法求解信号在不同载频、低频配置方案下各频率分量幅值，并计算相应的价值函数值。选取价值函数值最大者对应的载频、低频值作为网格搜索结果；
D.在网格搜索结果附近，采用二维二分搜索算法搜寻价值函数峰值位置，获得载频、低频精确估计值及对应的各频率分量幅值；
E.汇总估计结果，通过条件判断，判断接收信号是否为ZPW-2000信号以及信号是否受到干扰，输出载频、低频及提示信息。
3.根据权利要求2所述的一种基于非线性最小二乘法的ZPW-2000信号解调方法，其特征在于，所述步骤A的过程为：
A.1首先确认接收信号序列的长度N和采样频率fs；定义载频序列以及低频序列依次枚举载频和低频，形成144组编号为(p,q)的载频、低频配置方案分别构造对应的矩阵，·(p,q)表示编号为(p,q)的载频、低频配置方案下的数据；矩阵的参数配置方法参考式(3)和式(5)，其中的第一行元素采用CORDIC算法生成，的剩余元素由三角函数和角公式求得；由可得将144组与以全局变量的形式存入内存，供后续步骤(B.1)使用；
A.2定义载频序列 以及低频序列依次枚举载频和低频，形成144组编号为(p,q)的载频、低频配置方案对于编号(p,q)的配置方案，载频和低频的网格搜索间隔均为ΔωAII，定义如下：
式中，M(p,q)及其定义式中的NΔ和Δωd定义分别如下：
式中，为向下取整符号，意为四舍五入取整；以单位为Hz的实际频率值来描述，网格搜索中载频搜索频率值定义在区间内，低频搜索频率值定义在区间内；转化为数字频率，载频搜索频率值的定义如下：
其中，定义为：
当载频的网格搜索频率值确定后，即可以j为系数设定低频搜索值；j以及低频搜索值确定后各频率分量频率值的定义如下所示：
式中，·(p,q,i,j)表示编号(p,q)的载频、低频配置方案下编号(i,j)的网格搜索值对应的数据；和的定义如下：
将式(16)至式(20)代入式(5)，可得矩阵计算出因此，144组载频、低频配置方案中，每组方案对应一组M(p,q)，同时还包括组按次序以全局变量的形式将这些参数存入内存，供后续步骤(C.1)使用。
4.根据权利要求2所述的一种基于非线性最小二乘法的ZPW-2000信号解调方法，其特征在于，所述步骤B的过程为：
(B.1)令JBImax＝0，f0c＝-3000Hz，f1c＝-3000Hz；确认长度为N的接收信号序列x＝[x(0)x(1)…x(N-1)]T；依次提取内存中(A.1)存入的144组载频、低频配置方案对应的数据；其中，对于载频、低频配置方案其包含和采用高斯消元法求解方程：
将解得的αBI代入下式：
式中，Al，l＝2,3,4为各频率分量幅值，参考式(3)和式(8)由αBI求得；Athreshold为ZPW-2000信号包络幅值最小值，可实测获得；Ae为高信噪比下各频率分量幅值估计的最大误差，可实测获得，依据仿真结果，令Ae＝0.1；m为调制指数，其定义为r01为载频分量幅值与一次边频分量幅值比值的绝对值，定义如下：
若式(24)结果为真，则将αBI代入式(26)，得到该载频、低频配置方案的价值函数值：
否则令JBI＝0；此时，若JBI>JBImax，则更新粗估计价值函数及载频、低频粗估计值，即JBImax＝JBI，当完成144组载频、低频配置方案的枚举后，若f0c＝-3000(Hz)或f1c＝-3000Hz，则令并转至步骤E，否则转至步骤C。
5.根据权利要求2所述的一种基于非线性最小二乘法的ZPW-2000信号解调方法，其特征在于，所述步骤C的过程为：
C.1提取内存中A.2步骤存入的数据中满足条件的载频、低频配置方案对应的M(p,q)，计算：
ωshift＝ω0c-k0ΔωCI (27)
式中k0及其定义式中ω0c和ΔωCI定义分别如下：
定义向量xshift：
将序列x·xshift补零至M点，执行M点FFT得XM＝[XM(0) XM(1) XM(2) … XM(M-1)]；其中，·为向量点积符号；令JCImax＝0；依次提取A.2步骤存入数据中载频、低频配置方案组以(i,j)为编号的网格搜索数据；对于其中某一网格搜索编号(i,j)，计算η＝[η1 η2 η3 ... η2L]T，其中：
式中，Re(·),Im(·)分别为取实部符号和取虚部符号；ξl，l＝1,2,…,L，L＝5的定义如下：
ξ3＝XM(k0+i) (35)
采用高斯消元法计算方程：
将αCI代入下式：
JCI＝ηTαCI (39)
此时，若JCI>JCImax，则JCImax＝JCI，ω0g＝ω0c+iΔωCI，Δωg＝ΔωCI；当完成(2Lh+1)2组网格搜索值的枚举，转至步骤D。
6.根据权利要求2所述的一种基于非线性最小二乘法的ZPW-2000信号解调方法，其特征在于，所述步骤D的过程为：
D.1精确搜索由二维二分搜索算法实现。定义载频频率二分搜索初始区间的左边界ω0l、右边界ω0r及区间半长Δω0定义为：ω0l＝ω0g-Δωg，ω0r＝ω0g+Δωg，Δω0＝Δωg；搜索终止条件为2Δω0≤ωrp，ωrp为解调方法的精度要求；载频频率二分搜索步骤如下：
Step 1：以载频频率ω0l、ω0r为参数进行低频频率二分搜索，分别得到对应的价值函数值J0l、J0r。令ω0m＝(ω0l+ω0r)/2；
Step 2：以载频频率ω0m为参数进行低频频率二分搜索，得到对应的低频估计值ω0m_1以及对应的频率分量系数α0m和价值函数值J0m；令Δω0＝Δω0/2，ω0m_last＝ω0m；若J0l<J0r，则J0l＝J0m，ω0m＝ω0m+Δω0，否则J0r＝J0m，ω0m＝ω0m-Δω0；
Step 3：若2Δω0≤ωrp，则结束二分搜索，进而载频、低频及各频率分量系数的精确估计值分别为否则回到step 2；
低频频率二分搜索步骤与载频频率二分搜索基本一致；相对应地，低频频率二分搜索初始区间的左边界ω1l、右边界ω1r及区间半长Δω1定义为：ω1l＝ω1g-Δωg，ω1r＝ω1g+Δωg，Δω1＝Δωg；搜索终止条件为2Δω1≤ωrp；令代入的载频频率参数值为ω0p。低频频率二分搜索步骤如下：
Step 1：以ω0p为载频频率参数，ω1l、ω1r为低频频率参数，代入式(3)和式(5)，得到矩阵Z1l、Z1r，参考式(23)和式(26)，结合高斯消元法，得到对应的价值函数值J1l、J1r。令ω1m＝(ω1l+ω1r)/2；
Step 2：同样地，参考式(3)、式(5)、式(23)和式(26)求得载频频率ω0p，低频频率ω1m对应的价值函数值J1m以及频率分量系数α1m。令Δω1＝Δω1/2，ω1m_last＝ω1m；若J1l<J1r，则J1l＝J1m，ω1m＝ω1m+Δω1，否则J1r＝J1m，ω1m＝ω1m-Δω1；
Step 3：若2Δω1≤ωrp，则结束低频频率二分搜索，低频频率二分搜索输出的低频估计值，价值函数值以及频率分量系数分别为ω1m_last，J1m，α1m，否则回到step 2。
7.根据权利要求2所述的一种基于非线性最小二乘法的ZPW-2000信号解调方法，其特征在于，所述步骤E的过程为：
E.1若或则提示信号接收错误，结束解调过程，否则将代入下式：
式中，为各频率分量幅值，参考式(3)和式(8)由求得；为载频分量幅值与一次边频分量幅值的比值，参考式(24)和式(25)求得；若式(40)结果为真，则正常显示载频、低频信息结束解调过程，否则显示载频、低频信息的同时，提示信号受到干扰，结束解调过程。
CN201811451854.6A 2018-11-30 2018-11-30 一种基于非线性最小二乘法的zpw-2000信号解调方法 Expired - Fee Related CN109639612B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811451854.6A CN109639612B (zh) 2018-11-30 2018-11-30 一种基于非线性最小二乘法的zpw-2000信号解调方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811451854.6A CN109639612B (zh) 2018-11-30 2018-11-30 一种基于非线性最小二乘法的zpw-2000信号解调方法

Publications (2)

Publication Number Publication Date
CN109639612A true CN109639612A (zh) 2019-04-16
CN109639612B CN109639612B (zh) 2021-03-30

Family

ID=66070406

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811451854.6A Expired - Fee Related CN109639612B (zh) 2018-11-30 2018-11-30 一种基于非线性最小二乘法的zpw-2000信号解调方法

Country Status (1)

Country Link
CN (1) CN109639612B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114584901A (zh) * 2022-03-03 2022-06-03 西北工业大学 一种基于克罗内克分解的rls声反馈抑制算法

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1537770A (zh) * 2003-10-23 2004-10-20 北京交大思诺科技有限公司 轨道电路道碴泄漏电阻车载测试方法及其测试设备
CN102496064A (zh) * 2011-12-01 2012-06-13 北京交通大学 一种轨道不平顺的获取方法
CN102636693A (zh) * 2012-05-04 2012-08-15 重庆大学 一种结合fft与非线性最小二乘的谐波分析算法
CN102841381A (zh) * 2011-06-23 2012-12-26 中国石油天然气集团公司 一种基于分组线性拟合原理的单频干扰波压制方法
CN202814931U (zh) * 2012-09-22 2013-03-20 华南理工大学 一种基于频谱认知的自适应超声钢轨探伤装置
CN105634722A (zh) * 2015-12-28 2016-06-01 西安电子科技大学 一种mfsk伪装为跳频体制的抗截获方法
CN106524967A (zh) * 2016-11-07 2017-03-22 重庆理工大学 一种汽车轮心实际行驶位移测量与提取方法
CN106597408A (zh) * 2016-12-16 2017-04-26 重庆邮电大学 基于时频分析和瞬时频率曲线拟合的高阶pps信号参数估计方法
CN107329932A (zh) * 2017-05-08 2017-11-07 上海交通大学 基于非线性调频分量分解的时频域模态参数辨识方法
CN107612587A (zh) * 2017-06-20 2018-01-19 西安电子科技大学 一种用于跳频非合作通信中跳频信号的参数估计方法
CN108897917A (zh) * 2018-05-31 2018-11-27 南京东南建筑机电抗震研究院有限公司 一种用于高速铁路桥梁支座易损性评估的方法

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1537770A (zh) * 2003-10-23 2004-10-20 北京交大思诺科技有限公司 轨道电路道碴泄漏电阻车载测试方法及其测试设备
CN102841381A (zh) * 2011-06-23 2012-12-26 中国石油天然气集团公司 一种基于分组线性拟合原理的单频干扰波压制方法
CN102496064A (zh) * 2011-12-01 2012-06-13 北京交通大学 一种轨道不平顺的获取方法
CN102636693A (zh) * 2012-05-04 2012-08-15 重庆大学 一种结合fft与非线性最小二乘的谐波分析算法
CN202814931U (zh) * 2012-09-22 2013-03-20 华南理工大学 一种基于频谱认知的自适应超声钢轨探伤装置
CN105634722A (zh) * 2015-12-28 2016-06-01 西安电子科技大学 一种mfsk伪装为跳频体制的抗截获方法
CN106524967A (zh) * 2016-11-07 2017-03-22 重庆理工大学 一种汽车轮心实际行驶位移测量与提取方法
CN106597408A (zh) * 2016-12-16 2017-04-26 重庆邮电大学 基于时频分析和瞬时频率曲线拟合的高阶pps信号参数估计方法
CN107329932A (zh) * 2017-05-08 2017-11-07 上海交通大学 基于非线性调频分量分解的时频域模态参数辨识方法
CN107612587A (zh) * 2017-06-20 2018-01-19 西安电子科技大学 一种用于跳频非合作通信中跳频信号的参数估计方法
CN108897917A (zh) * 2018-05-31 2018-11-27 南京东南建筑机电抗震研究院有限公司 一种用于高速铁路桥梁支座易损性评估的方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
戈壁舟等: "ZPW-2000移频信号二分搜索频域解调法及窗函数的选择", 《兰州交通大学学报》 *
秦晓荣: "浅析ZPW2000A轨道电路工作频率的干扰", 《信息通信》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114584901A (zh) * 2022-03-03 2022-06-03 西北工业大学 一种基于克罗内克分解的rls声反馈抑制算法

Also Published As

Publication number Publication date
CN109639612B (zh) 2021-03-30

Similar Documents

Publication Publication Date Title
CN109167746B (zh) 连续波与脉冲信号快速识别装置
US9806914B1 (en) Successive signal interference mitigation
Zeng et al. Automatic modulation classification of radar signals using the Rihaczek distribution and Hough transform
CN111693944B (zh) 雷达有源干扰信号参数提取与干扰样式识别方法及装置
CN107290589A (zh) 基于短时分数阶傅里叶变换的非线性信号时频分析方法
CN108548957B (zh) 基于循环调制频谱和分段互相关相结合的双谱分析方法
CN107102255A (zh) 单一adc采集通道动态特性测试方法
CN110133632B (zh) 一种基于cwd时频分析的复合调制信号识别方法
CN109839618B (zh) 低信噪比雷达信号识别方法、计算机可读存储介质及系统
CN106357574A (zh) 基于顺序统计量的bpsk/qpsk信号调制盲识别方法
CN105680905B (zh) 一种适用于任意调制度的fm、pm信号载波捕获方法
CN107171994B (zh) 无线电引信信号识别和重构系统及方法
CN111474524A (zh) 一种雷达干扰装备干扰效果监测与决策支持系统
CN107561497B (zh) Fsk及多种非线性调频信号识别及参数估算方法
CN109639612B (zh) 一种基于非线性最小二乘法的zpw-2000信号解调方法
CN110187313A (zh) 基于分数阶Fourier变换的雷达信号分选识别方法及装置
CN109342813B (zh) 一种基于dft和二分法的正弦信号频率估计方法
CN104363194A (zh) 基于波形变换的psk调制识别方法
CN108572352B (zh) 一种基于欠采样的相位编码信号的参数估计方法
CN109379310A (zh) 一种基于Rife-Quinn综合的MPSK信号载频估计方法
CN107835036A (zh) 非合作跳频信号破解方法
CN109067680A (zh) 一种基带信号的载波频偏估计方法及其装置
CN101588191B (zh) 无线电信号认知方法及设备
CN103905360B (zh) 加入“判极”操作的非合作bpsk信号解码方法
CN109682492B (zh) 基于频域高斯拟合的频率估计方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210330

Termination date: 20211130