CN112769725B - 基于全相位频谱纠正的Costas序列时频联合同步方法 - Google Patents

Abstract

本发明涉及数字信号处理领域，具体涉及一种基于全相位频谱纠正的Costas序列时频联合同步方法，该方法对使用短时傅里叶变换和滑动相关的离散频谱纠正算法所存在的问题进行改进，改进主要包括：本发明提出了一种基于迭代寻优的解决方案，当实际的频偏是频谱分辨率的整数倍时会出现较大的误差，本发明使用估计的误差较大的估计值对信号进行频偏纠正与时延纠正，然后再次进行估计，通过迭代的方式使得信号的频偏不是特殊值。针对问题二，由于频谱泄露引起的误差，本发明使用全相位傅里叶变化代替短时傅里叶变换中的FFT来最大可能的减小频谱泄露现象。通过改进以上两点来提高原有算法的估计精度。

Description

基于全相位频谱纠正的Costas序列时频联合同步方法

技术领域

本发明涉及数字信号处理领域，具体涉及一种基于全相位频谱纠正的Costas序列时频联合同步方法。

背景技术

短波MFSK通信被广泛应用在低速率最低保障的应用场景中。对于通信系统而言，接收信号的前提是检测信号以及实现收发端同步，同步性能的好坏直接影响通信系统的性能。出现同步误差或失去同步会导致通信系统性能下降或通信中断。

在短波MFSK通信系统中，常用的信号同步方法是数据辅助类同步算法。数据辅助类同步算法，通过在数据序列中插入训练序列完成，接收端通过本地序列寻找数据序列中的训练序列来进行同步。常用的训练序列有自相关性良好的伪随机序列等。常用的同步算法有滑动相关法等。由于Costas序列在二维矩阵表示上具有优秀的二维自相关特性，因此基于Costas序列的短波MFSK通信时频联合同步算法使用Costas序列作为同步序列，通过短时傅里叶变换讲时域信号处理成时频谱的形式，使用Costas序列作为本地序列进行二维滑动相关。得益于Costas序列二维矩阵的优秀自相关特性，该同步算法具有优异的抗噪声性能和抗频偏性能。基于Costas序列的短波MFSK通信时频联合同步算法频率同步精度受限于频谱分辨率，时间同步精度受限于码元时间长度，因此可以通过过采样来提高频率估计精度，通过改变短时傅里叶变换中窗函数的步进长度来提高时间同步精度。除此之外，可以使用离散频谱分析算法提高估计精度。但由于离散频谱分析算法在使用中的两点问题，造成了基于Costas序列的短波MFSK通信时频联合同步算法的估计精度没有达到最优。问题描述如下：

(1)在实际的频率偏移值f是频谱分辨率的整数倍时(可以为0，即无频偏)，会出现一个较大的误差。当无频偏时，在无噪声干扰的情况下，假设经过DFT后的峰值谱线其谱线位置应该为x处，这种情况下，x处的谱线也就是峰值处即为正确的频率位置，且峰值谱线的且x-1处和x+1处的值为0，重心法则可以估计出正确的频偏。但由于受到噪声的干扰，会造成x-1和x+1处的值不为0，使离散频谱分析算法造成较大误差。

(2)在第一点的分析中没有考虑时延，下面对时延造成的影响进行分析。如果信号存在时延，则其采样后的信号会存在一段只含有噪声的信号，这种情况下，使用长度为单个码元时间的窗函数进行短时傅里叶变换会导致当前码元的部分信息存在于下一码元时间，在MFSK信号中，当两个码元的频率间隔较小时，由于频谱泄露的影响，就会对离散频谱分析算法所使用的主瓣谱线造成严重的干扰，造成较大误差。

综上所述，虽然基于Costas序列的短波MFSK通信时频联合同步算法具有较高的估计精度，但仍存在可以改进的空间，对频谱纠正算法的改进，将进一步提高该算法的同步精度。

发明内容

为了解决上述问题，本发明提供一种基于全相位频谱纠正的Costas序列时频联合同步方法。

一种基于全相位频谱纠正的Costas序列时频联合同步方法，包括以下步骤：

S1、接收端收到时域信号x(N)后，按照码元时间将时域信号x(N)进行分割，得到分割信号为x₁(n)、x₂(n)…x_m(n)，分割信号用x_i(n)表示，其中i＝1,2,……,m，m表示分割信号的数量；

S2、将x_i(n)进行翻转变换后加窗进行全相位傅里叶变换apFFT，完成全相位傅里叶变换过程后得到分割信号x_i(n)的全相位傅里叶变换X₁(W)--X_m(W)，进行模值计算得到分割信号x_i(n)的频谱|X₁(W)|²--|X_m(W)|²，用|X_i(W)|²表示；

S3、将分割信号的频谱|X₁(W)|²--|X_m(W)|²合并成一个时频谱，并将合并后的时频谱表示为时频矩阵；

S4、使用滑动窗口法将本地Costas序列在时频矩阵上进行滑动自相关，得到二维自相关函数，通过自相关函数的峰值坐标获得粗同步结果(即时延粗估计值和频偏粗估计值)，根据粗估计值对全相位傅里叶变换后的离散频谱进行纠正，得到时延估计值和频偏估计值；

S5、判断：判断步骤S4中二维自相关函数的峰值幅度与上次迭代相比是否减小或者迭代次数是否达到最大限制；如果二维自相关函数的峰值幅度减小或者迭代次数达到最大限制，则输出上次迭代中步骤S4的时延估计值和频偏估计值作为算法输出结果；如果没有减小，则迭代次数加一，执行步骤S6；

S6、将步骤S4的时延估计值和频偏估计值作为纠正值对步骤S1的时域信号x(N)进行时延和频偏纠正，获得纠正后的时域信号x(N)，返回执行步骤S1。

进一步的，步骤S4具体包括：

S41、使用本地Costas序列在步骤S3所获得时频矩阵上进行二维自相关，获得二维自相关函数，取二维自相关函数的峰值坐标的横坐标作为时延粗估计值，取峰值坐标的纵坐标作为频偏粗估计值；

S42、根据时延粗估计值，取时延粗估计值对应的码元(时延粗估计值确定)时间内的频谱进行离散频谱分析，获得该频谱峰值坐标；通过重心法则和所使用的窗函数主瓣函数表达式进行频率粗估计，获得时延估计值和频偏估计值。

进一步的，步骤S6中，纠正公式包括：

x(n)=x(n+t_es)·e^-j2πΔft

其中，x(n)为纠正后的时域信号，t_es为时延粗估计值，Δf为频偏估计值。

本发明的有益效果：

1.本发明方法使用全相位傅里叶变化代替传统方法中所使用的短时傅里叶变换STFT，通过apFFT能够最大可能的减小MFSK信号中由于时延所造成的相邻码元频谱泄漏和互相干扰问题，从而可以有效避免后续由于频谱泄漏和互相干扰造成的离散频谱分析误差较大的问题。

2.传统的基于Costas序列的时频联合同步算法是基于短时傅里叶变换，通常通过过采样提高频谱分辨率，但这种方式会导致公式中N的数值增大，使得算法的抗噪声能力下降；此外，在实际的频率偏移值f是频谱分辨率的整数倍时，会出现较大误差，本发明通过迭代对该误差进行频偏纠正与时延纠正。

附图说明

下面结合附图和具体实施方式对本发明做进一步详细的说明。

图1为本实施例中的一种基于apFFT的时频联合同步算法流程图；

图2为本实施例中的一种基于迭代寻优的离散频谱纠正算法流程图；

图3为频谱泄漏影响说明图；

图4为本实施例中的一种离散频谱纠正示意图1；

图5为本实施例中的另一种离散频谱纠正示意图2。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本实施例针对现有的基于Costas序列的短波MFSK通信时频联合同步算法中存在的估计精度问题，对方法进行改进，实现了离散频谱纠正算法的精度提升。改进之处主要包括使用apFFT对基于Costas序列的短波MFSK通信时频联合同步算法的改进和基于迭代寻优的离散频谱纠正算法改进。

实施例一

本实施例提供了一种基于全相位频谱纠正的Costas序列时频联合同步方法，如图1-2所示，在一种优选实施方式中，其包括以下步骤：

步骤S1，接收端收到时域信号x(N)后，按照码元时间将时域信号x(N)进行信号分割，将时域信号x(N)变为m个表示码元信息的时域信号，得到分割信号为x₁(n)、x₂(n)…x_m(n)，简称为分割信号x_i(n)，其中i＝1,2,……,m，m表示分割信号的数量。

步骤S2，将x_i(n)进行翻转变换后加窗进行全相位傅里叶变换，完成apFFT过程后得到分割信号x_i(n)的全相位傅里叶变换X₁(W)--X_m(W)，进行模值计算得到分割信号x_i(n)的频谱|X₁(W)|²--|X_m(W)|²，用|X_i(W)|²表示。

具体地，步骤S2包括：将分割信号x_i(n)复制并翻转后再与信号x_i(n)组合，得到信号x_i(2n-1)。

然后使用卷积窗w对组成信号x_i(2n-1)进行加权，将间隔为N的数据两两进行重叠相加，变为信号y_i(n)，N代表的是数字信号整体整体的长度。

对信号y_i(n)进行加窗DFT，完成全相位傅里叶变换运算，得到分割信号x_i(n)的全相位傅里叶变换X₁(W)--X_m(W)。随后对加窗DFT后的y_i(n)信号进行模值计算，信号y_i(n)经过模值计算后，得到分割信号x₁(n)--x_m(n)的频谱|X₁(W)|²--|X_m(W)|²，信号x_i(n)的频谱用|X_i(W)|²表示。

本步骤在使用码元时间长度分割时域信号后，对分割后的多个时域信号单独进行全相位傅里叶变化，得到多个时域信号的频谱，然后再使用将多个时域信号的频谱合并的方式得出信号x(N)的时频谱。该步骤的有益效果包括：为了减小由于频谱泄露的影响，使用全相位傅里叶变化代替传统短时傅里叶变换(STFT)中的FFT过程，通过apFFT能够最大可能的减小MFSK信号中由于时延所造成的相邻码元频谱泄漏和互相干扰问题，从而可以有效避免后续由于频谱泄漏和互相干扰造成的离散频谱分析误差较大的问题。

步骤S3，将分割信号x_i(n)的频谱|X₁(W)|²--|X_m(W)|²合并成一个时频谱，得到合并后的时频谱，并用时频矩阵表示合并后的时频谱。

步骤S4，使用滑动窗口法将本地Costas序列在时频矩阵上进行滑动自相关，通过自相关函数的峰值坐标获得粗同步结果。具体地，使用本地Costas序列在步骤S3所获得时频矩阵上进行二维自相关，获得二维自相关函数后，取二维自相关函数的峰值坐标的横坐标作为时延粗估计值t_es，取二维自相关函数的峰值坐标的纵坐标作为频偏粗估计值。

由于频偏f的影响，x为经过apFFT变换的频谱峰值位置，如图4所示，x₀为真实谱峰位置，x为经过apFFT变换的频谱峰值位置，因此，上述直接取二维自相关函数峰值坐标的纵坐标作为频偏粗估计值存在一定误差，因此需要对频偏进行进一步纠正，得到对频偏粗估计值纠正后的频偏估计值。具体地，根据时延粗估计值，通过t_es取相应码元时间内的频谱进行离散频谱分析，获得该频谱峰值坐标，通过重心法则和所使用的窗函数主瓣函数表达式进行频率估计，得到频偏估计值Δf。在这一步中获得时延估计值t_es和频偏估计值Δf。

步骤S5，判断：通过判断步骤S4中二维自相关函数的峰值幅度与上次迭代(如果为初次，则上次迭代的二维自相关函数的峰值幅度默认为0)相比是否减小，或者迭代次数是否达到最大限制。如果减小或者迭代次数达到最大限制，则输出上次迭代中步骤S4的时延估计值和频偏估计值作为算法输出结果；如果没有减小，则迭代次数加一，执行步骤S6。

具体地，在步骤S1中得到二维自相关函数，通过比较此次的二维自相关函数的幅度Max值,定义为Max_i，i代表迭代次数。将Max_i与Max_i-1(如果是第一次执行步骤S2,Max_i-1默认为0，即Max₁默认为0)相比，如果二维自相关函数的峰值幅度与上次迭代相比减小，即Max_i＞Max_i-1，则说明此次的峰值谱线位置与上次相比更加靠近真实峰值，其峰值能量与上次相比变大，代表频偏纠正为正效果，所以则迭代次数加一，继续执行后续步骤S6。如果二维自相关函数的峰值幅度与上次迭代相比没有减小，即Max_i＜＝Max_i-1，则说明此次的峰值谱线位置与上次相比更加远离真实峰值，其峰值能量与上次相比变小，代表频偏纠正为负效果，应该及时跳出，防止下次迭代出现纠正为正效果的情况，出现无限迭代的情况，直到最大迭代次数。时延τ-t_es对Max_i的影响与频偏f-Δf相同，所以Max_i＜＝Max_i-1时，跳出迭代，并输出上次迭代步骤S4的时延估计t_es和频偏估计Δf作为最终估计结果。

在本实施例中通过前后两次迭代的自相关函数峰值作为是否结束迭代的依据。通过比较此次的Max值,定义为Max_i，i代表迭代次数。将Max_i与Max_i-1相比，如果Max_i＞Max_i-1，则继续进行迭代提高算法精度。如果Max_i＜＝Max_i-1则说明此次的峰值谱线位置与上次相比更加远离真实峰值，其峰值能量与上次相比变小，代表频偏纠正为负效果，应停止迭代，输出结果。算法以倒数第二次迭代作为输出结果，保证算法的精度。

步骤S6，将步骤S4的时延估计值t_es和频偏估计值Δf作为纠正值对步骤S1的时域信号x(N)进行时延和频偏纠正，获得纠正后的时域信号x(N)，对纠正后的时域信号x(N)执行步骤S1-S5。

本发明方法以迭代的方式来解决离散频谱分析中，频偏为频谱分辨率的整数倍时所造成的的大误差，每次迭代都按照估计值(不管实际精度如何)对原始信号进行纠正，通过不断地纠正，使得频偏处于特殊值的概率越来越小，从而使得离散频谱分析算法的估计精度达到最优。

综上，本发明的一种基于全相位频谱纠正的Costas序列时频联合同步方法，能够有效改善时频联合同步算法同步精度呈现周期性的问题，且能够进一步提高算法的抗噪声性能。

实施例二

本实施例中假设码元时间为T，采样率为fs，时域信号频偏为f，时延为τ。

一种基于全相位频谱纠正的Costas序列时频联合同步方法，包括但不限于如下步骤：

(1)、将时域信号按照码元时间T进行分割，分割信号为x₁(n)、x₂(n)…x_m(n)。由于时延的影响，前一个码元的部分信息会存在与后一个分割信号内。以x₂(n)为例，如果没有时延，则x₂(n)所表示的是第二个码元符号的时域波形，此时x₂(n)为：

x₂(n)＝A₂ sin[2πf₂(n-τ)] (1)

其中，A₁表示第一个码元时间内的信号幅度，A₂表示第二个码元时间内的信号幅度，n为时间离散表示，τ表示时延，f₁为码元1所代表的频率，f₂为码元2所代表的频率此时对x₂(n)进行FFT。频率f₁自身所产生的频谱泄漏现象并不会对后续估计产生影响。当存在时延时，第一个码元符号的部分信息会延后到x₂(n)中，即

x₂(n)=A₁ sin[2πf₁(n-τ)]+A₂ sin[2πf₂(n-τ)] (2)

此时对x₂(n)进行傅里叶变换FFT后，频率f₁所造成的频谱泄漏现象会极大的干扰频率f₁的主瓣谱线，如图3所示。此时的干扰会严重的影响后续的离散频谱纠正。

为了尽可能减小由于上述频谱泄漏现象造成的频谱纠正问题，本实施例在FFT之前对信号进行处理(主要包括信号转换、窗函数加权等处理，详见下述描述)，并使用全相位傅里叶变换apFFT代替传统的FFT，尽可能的减小这种情况下频谱泄漏造成的影响。具体描述如下：

(2)、将分割信号进行信号反转后与原信号拼接，组合成新的信号x_i(2n-1)，使用卷积窗w_c对信号x_i(2n-1)进行加权，将间隔为N的数据两两进行重叠相加，N在数值上为分割后信号x_i(n)的数据点数。w_c为长为n的前窗g₁(n)和翻转的后窗g₂(n)的卷积。

w_c(n)=g₁(n)*g₂(-n) (3)

随后将间隔为N的数据两两进行重叠相加，并加窗函数R(n)变为信号y_i(n)，y_i(n)的表达式如下：

y_i(n)=[w_c(n)x_i(n)+w_c(n-N)x_i(n-N)]R(n) (4)

(3)、将(2)中的y_i(n)进行DFT，就完成了对分割信号x_i(n)的全相位傅里叶变换，然后对信号进行模值计算。定义Y_i(W)为信号x_i(n)的全相位傅里叶变换，则x_i(n)的谱图为Y_i(W)的模值，表达式如下：

x_i(n)→|Y_i(W)|² (5)

箭头代表|X_i(W)|²|Y_i(W)|²为时域信号x_i(n)的频域映射。

将分割信号x_i(n)的频谱|Y₁(W)|²--|Y_m(W)|²合并成一个时频谱，得到合并后的时频谱，并用时频矩阵表示合并后的时频谱，则时频信号x(n)的时频谱表示如下：

上式中，双箭头代表信号x(n)的时频域映射，为矩阵形式。

然后使用本地Costas序列进行滑动同步，获得二维自相关函数的峰值Max，峰值坐标为(t_es,f_es)，t_es为时延粗估计值，f_es为频率粗估计值。

(4)时延粗估计值t_es的分辨率与(1)中的码元时间T相同，所以根据t_es的值可以获得x(n)的时频谱中某个码元符号的频谱|Y_i(W)|²，f_es所表示的就是该码元符号的频谱|Y_i(W)|²的峰值处，其谱图如图4所示。由于频偏f的影响，x为经过apFFT变换的频谱峰值位置，x₀为真实谱峰位置。如果图1中的窗函数为矩形窗，则可以定义其功率谱函数Y_i(x)表达式如下：

根据重心法则，在n→∞时：

其中，n为频谱主瓣内谱线数的一半，根据公式(7)(8)可以计算出频谱函数真实谱峰位置x₀的位置坐标为：

设采样频率为f_s，FFT点数为N，定义k为谱线峰值坐标，则其频偏估计值为：

可以看出，所估计的Δf精度受到f_s/N的限制，也就是受到频谱分辨率的限制。传统的基于Costas序列的时频联合同步算法是基于短时傅里叶变换，通常通过过采样提高频谱分辨率，但过采样会导致公式中N的数值增大，使得算法的抗噪声能力下降。且在实际的频率偏移值f是频谱分辨率的整数倍时(可以为0，即无频偏)，会出现一个较大的误差。如图5所示，当无频偏时，在无噪声干扰的情况下，x处的谱线也就是峰值处即为正确的频率位置，且x-1处和x+1处的值为0，根据公式(10)可以估计出正确的频偏。但由于受到噪声的干扰，会造成x-1和x+1处的值不为0，再使用公式(10)进行计算就会造成较大误差。

(5)针对上述(4)中可能存在的误差，通过图2所示流程解决。当接收端接收到时域信号后，开始经过图2流程处理，如果频偏如(4)中所描述，是频谱分辨率的整数倍时(可以为0，即无频偏)，经过公式(6)会得出一个误差可能较大的频偏值Δf，因此，本实施例对信号进行频偏纠正与时延纠正(时延纠正的值为t_es)，纠正公式为：

x(n)=x(n+t_es)·e^-j2πΔft (11)

其中，x(n)为图1中的时域信号。则经过公式(11)后，时域信号频偏为f-Δf，时延为τ-t_es。频偏f-Δf依然处于图5所示的特殊位置处的概率将大大减小。然后再次从步骤(1)进行，在进行到步骤(3)时，获得此次的Max值，定义为Max_i，i代表迭代次数且满足i>1。将Max_i与Max_i-1相比，通过以下方案判断。

其判断依据是，如果Max_i＞Max_i-1，则说明此次的峰值谱线位置与上次相比更加靠近真实峰值，其峰值能量与上次相比变大，代表频偏纠正为正效果，所以继续执行后续步骤到(5)。如果Max_i＜＝Max_i-1则说明此次的峰值谱线位置与上次相比更加远离真实峰值，其峰值能量与上次相比变小，代表频偏纠正为负效果，应该及时跳出，防止下次迭代出现纠正为正效果的情况，出现无限迭代的情况，直到最大迭代次数。时延τ-t_es对Max_i的影响与频偏f-Δf相同，所以Max_i＜＝Max_i-1时，跳出迭代，并输出上次迭代步骤(3)的t_es和步骤(4)的Δf作为最终估计结果。

当介绍本申请的各种实施例时，词语“包括”、“包含”和“具有”都是包括性的并意味着除了列出的步骤之外，还可以有其它步骤。

需要说明的是，本领域普通技术人员可以理解实现上述方法实施例中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法实施例的流程。其中，所述存储介质可为磁碟、光盘、只读存储记忆体(Read-0nly Memory，ROM)或随机存储记忆体(RandomAccess Memory，RAM)等。

以上所述仅是本申请的具体实施方式，应当指出，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (3)

1.一种基于全相位频谱纠正的Costas序列时频联合同步方法，其特征在于，包括以下步骤：

S1、接收端收到时域信号x(N)后，按照码元时间将时域信号x(N)进行分割，得到分割信号为x₁(n)、x₂(n)…x_m(n)，分割信号用x_i(n)表示，其中i＝1,2,……,m，m表示分割信号的数量；

S2、将x_i(n)进行翻转变换后加窗进行全相位傅里叶变换apFFT，完成全相位傅里叶变换过程后得到分割信号x_i(n)的全相位傅里叶变换X₁(W)--X_m(W)，进行模值计算得到分割信号x_i(n)的频谱|X₁(W)|²--|X_m(W)|²，用|X_i(W)|²表示；所述将x_i(n)进行翻转变换包括将分割信号x_i(n)复制并翻转后再与信号x_i(n)组合；

S3、将分割信号的频谱|X₁(W)|²--|X_m(W)|²合并成一个时频谱，并将合并后的时频谱表示为时频矩阵；

S4、使用滑动窗口法将本地Costas序列在时频矩阵上进行滑动自相关，得到二维自相关函数，通过自相关函数的峰值坐标获得粗同步结果，根据粗估计值对全相位傅里叶变换后的离散频谱进行纠正，得到时延估计值和频偏估计值；

S5、判断：判断步骤S4中二维自相关函数的峰值幅度与上次迭代相比是否减小或者迭代次数是否达到最大限制；如果二维自相关函数的峰值幅度减小或者迭代次数达到最大限制，则输出上次迭代中步骤S4的时延估计值和频偏估计值作为算法输出结果；如果没有减小，则迭代次数加一，执行步骤S6；

S6、将步骤S4的时延估计值和频偏估计值作为纠正值对步骤S1的时域信号x(N)进行时延和频偏纠正，获得纠正后的时域信号x(N)，返回执行步骤S1。

2.根据权利要求1所述的一种基于全相位频谱纠正的Costas序列时频联合同步方法，其特征在于，步骤S4具体包括：

S41、使用本地Costas序列在步骤S3所获得时频矩阵上进行二维自相关，获得二维自相关函数，取二维自相关函数的峰值坐标的横坐标作为时延粗估计值，取峰值坐标的纵坐标作为频偏粗估计值；

S42、根据时延粗估计值，取时延粗估计值对应的码元时间内的频谱进行离散频谱分析，获得该频谱峰值坐标；通过重心法则和所使用的窗函数主瓣函数表达式进行频率粗估计，获得时延估计值和频偏估计值。

3.根据权利要求1所述的一种基于全相位频谱纠正的Costas序列时频联合同步方法，其特征在于，步骤S6中，纠正公式包括：

x(n)＝x(n+t_es)·e^-j2πΔft

其中，x(n)为纠正后的时域信号，t_es为时延粗估计值，Δf为频偏估计值。

Images