技术领域
[0001] 本
发明涉及立方相位信号参数估计技术领域，更具体地，本发明涉及雷达目标检测和参数估计技术领域中的一种基于循环移位的立方相位信号参数估计方法。
背景技术
[0002] 在雷达、声纳、无线通信、地质勘探和医学成像等实际工程技术应用中，多数信号具有非平稳特性，其相位也具有连续瞬时性。因此，多项式相位信号成为该类领域处理信号最为常见和广泛使用的一种信号模型。
[0003] 多项式相位信号的相位参数包含着重要信息。在雷达
信号处理中，一次相位、二次和三次相位参数分别对应着高速目标的初始速度、
加速度及加速度变化率，其估计的准确性对目标的检测和成像具有至关重要的作用。在无线通信中，由于传播介质存在物理扰动，接收系统与目标之间的相对运动产生的
多普勒效应会使信号的
频率发生变化。在
生物界及生活中的信号，如蝙蝠信号、脑电波信号等都可近似视为多项式相位信号，
地震波探测和医学成像等其他领域中多项式相位信号也具有应用价值。因此，多项式相位信号已成为非平稳信号处理领域中的重要研究对象，其研究重点在于基于立方相位信号的检测和参数估计。说明立方相位信号的定义和概念。
[0004] 目前，基于时频分析方法的立方相位信号参数估计方法大致有离散调频率傅里叶变换(Discrete Chirp Fourier Transform，DCFT)、多项式Winger-Ville分布(Polynomial Winger-Ville Distribution，PWVD)、高阶模糊函数法(High-Order Ambiguity Function，HAF)以及基于HAF的乘积型高阶模糊函数(Product High-Order Ambiguity Function，PHAF)等
算法。然而，PWVD算法和HAF算法需要利用高阶多重变换求取信号的高阶瞬时自相关函数，再进行傅里叶变换，先估计出最高阶系数，再通过多次解调频(Dechirp)降
低信号相位阶数，依次估计其他相位系数，如此重复，直到估计出所有的相位参数。该类方法的
缺陷在于对于多分量的立方相位信号参数是从高阶到低阶是顺序估计的，低阶相位参数的估计性能依赖于高阶参数估计的精确度，这会引起严重的误差传递(Error Propagation)效应。另外，此类方法都是对非均匀
采样数据采用傅里叶变换，无法利用FFT快速的完成相应的傅里叶变换，而且会对参数的估计
精度产生一定的影响。针对立方相位信号的DCFT算法不会产生交叉项问题，其不足之处在于它是通过二维搜索实现调频率和二次调频率的估计，搜索过程将导致庞大的运算量，增加系统的复杂度。
发明内容
[0005] 本发明的目的在于针对上述
现有技术的不足，提出基于循环移位的立方相位信号参数估计方法。该方法弥补了传统多项式相位信号参数估计算法中多个参数依次估计、搜索过程复杂、计算量大等缺陷，利用非均匀快速傅里叶变换(Non-uniform Fast Fourier Transform，NUFFT)实现信号的时频分布以及循环移位的方法，实现了多个参数同时估计，同时利用NUFFT算法和循环移位操作大大减小了检测过程中的计算量。
[0006] 实现本发明的思路是：首先给出离散立方相位信号的参数表达式，然后利用非均匀快速傅里叶变换求取该立方相位信号的时频(Time-Doppler Frequency)分布，再对已得到的时频分布转换其参数空间并进行循环移位累加操作，最后利用峰值检测同时完成调频率和二次调频率的联合估计。
[0007] 为实现上述技术目的，本发明采用如下技术方案予以实现。
[0008] 基于循环移位的立方相位信号参数估计方法包括以下步骤：
[0009] 步骤1，获取单分量立方相位信号z(n)，n表示离散时间变量；
[0010] 步骤2，定义所述单分量立方相位信号的瞬时自相关函数Rz(n,m)，Rz(n,m)＝z(n+m)z(n-m)，其中，m表示离散延时变量；
[0011] 步骤3，对所述单分量立方相位信号的瞬时自相关函数Rz(n,m)沿离散延时变量m的坐标轴进行非均匀快速傅里叶变换，得出所述单分量立方相位信号的时频分布Wz(σ,n)，得出所述单分量立方相位信号的时频分布的模值Dσ,n，σ表示频率变量；
[0012] 步骤4，建立参数空间n-σ-D，参数空间n-σ-D为大小为M×N的矩阵，N表示所述单分量立方相位信号的时域长度，M表示所述单分量立方相位信号的频域长度；参数空间n-σ-D第l行第k列的取值为 nk表示离散时间变量n的第k个取值，k＝1,2,...,N，σl表示离散频率变量σ的第l个取值，l＝1,2,...,M；
[0013] 通过对参数空间n-σ-D的每列数据进行循环移位，建立参数空间β-α-Ψ，参数空间β-α-Ψ是大小为M×N的矩阵；令i＝1,2,...,N，令j＝1,2,...,M；参数空间β-α-Ψ第j行第i列的取值 应为：满足方程σ＝αj+βin的所有二维坐标(n,σ)对应的时频分布的模值之和，αj表示截距α的第j个取值，βi表示斜率β的第i个取值；
[0014] 步骤5，得出参数空间β-α-Ψ元素最大值对应的行序号j(max)和参数空间β-α-Ψ元素最大值对应的列序号i(max)；得出所述单分量立方相位信号的调频率的估计值 以及所述单分量立方相位信号的二次调频率的估计值 为截距α的第j(max)个取值αj(max)，为斜率β的第i(max)个取值βi(max)。
[0015] 本发明的有益效果为：
[0016] 1)本发明利用非均匀快速傅里叶变换(NUFFT)获取信号的时频分布，实现信号
能量的有效积累的，可以有效减小计算量。
[0017] 2)本发明利用转换数据空间的方法，实现了立方相位信号的多个参数同时估计，避免了单个参数逐个估计过程中出现的误差传递，摆脱了复杂的二维谱峰搜索过程，克服了现有技术计算量大、过程复杂的不足。
附图说明
[0018] 图1为本发明的基于循环移位的立方相位信号参数估计方法的
流程图；
[0019] 图2a为本发明
实施例中直线检测的第一原理示意图；
[0020] 图2b为本发明实施例中直线检测的第一原理示意图；
[0021] 图3为仿真实验中得出的单分量立方相位信号的时频分布的示意图；
[0022] 图4为图3中所示的时频分布的等高线图；
[0023] 图5为仿真实验中得出新的参数空间的示意图；
[0024] 图6为仿真实验中得出的二次调频率-调频率分布示意图；
[0025] 图7为图6的等高线示意图。
具体实施方式
[0026] 下面结合附图对本发明作进一步说明：
[0027] 参照图1，为本发明的基于循环移位的立方相位信号参数估计方法的流程图。该基于循环移位的立方相位信号参数估计方法包括以下步骤：
[0028] 步骤1，获取单分量立方相位信号，得出所述单分量立方相位信号的离散信号表达式。
[0029] 其具体步骤为：
[0030] 在实际生活中获取信号，如雷达信号、脑电波信号、
地震波信号。将实际生活中获取的信号(实际信号)近似用单分量立方相位信号进行表示。
[0031] 本发明实施例中，单分量立方相位信号的连续信号表达式为：
[0032]
[0033] 其中，A为信号幅度，t为连续时间变量，f表示上述实际生活中获取的信号的中心频率(Centroid Frequency，CF)，α表示上述实际生活中获取的信号的调频率(Chirp Rate，CR)，β表示上述实际生活中获取的信号的二次调频率(Quadratic Chirp Rate，QCR)。
[0034] 根据对实际生活中获取的信号进行采样接收时的
采样频率Fs对上述单分量立方相位信号进行离散化，上述单分量立方相位信号的离散信号表达式为：
[0035]
[0036] 其中，n表示离散时间变量(取整数值)，N表示上述实际生活中获取的信号的长度。
[0037] 步骤2：定义该单分量立方相位信号的瞬时自相关函数Rz(n,m)，Rz(n,m)＝z(n+m)z(n-m)，其中，m表示延时变量。
[0038] 其具体步骤为：
[0039] 为了分析该单分量立方相位信号的时频特性，在此其定义该单分量立方相位信号的瞬时自相关函数Rz(n,m)，Rz(n,m)＝z(n+m)z(n-m)，其中，m表示延时变量，其长度为M。
[0040] 由述单分量立方相位信号的离散信号表达式可知，单分量立方相位信号的瞬时自相关函数Rz(n,m)为：
[0041]
[0042] 其中，m表示离散延时变量。可以看出，离散时间变量n和延时变量m存在线性耦合。这种耦合性正是该立方相位信号检测和参数估计的切入点。
[0043] 步骤3，对上述单分量立方相位信号的瞬时自相关函数Rz(n,m)沿离散延时变量m的坐标轴进行非均匀快速傅里叶变换(NUFFT)，得出上述单分量立方相位信号的时频分布(Time-Doppler Frequency分布)Wz(σ,n)、以及上述单分量立方相位信号的时频分布的模值Dσ,n，σ表示离散多普勒频率变量。
[0044] 其具体步骤为：
[0045] 将上述单分量立方相位信号的时间-多普勒频率分布Wz(σ,n)定义为：
[0046]
[0047] 其中，σ表示离散多普勒频率变量，m的取值范围为[0,(N-1)/2]。
[0048] 将 带入上式，计算得到
[0049]
[0050] 对其求模可得：
[0051] Dσ,n＝|Wz(n,σ)|＝|A2|*δ(σ-α-βn)
[0052] 其中，Dσ,n表示Wz(n,σ)的模值。
[0053] 从上式可以明显看出，时频分布Wz(n,σ)的模值Dσ,n可以看出沿着直线σ-α-βn＝0(在直线σ-α-βn＝0，σ是自变量n的函数)分布的冲击线谱，即时频分布的能量集中出现在表示信号瞬时频率变化律的直线上，该直线的截距和斜率分别正是上述单分量立方相位信号的调频率α和二次调频率β。
[0054] 步骤4，建立参数空间n-σ-D，参数空间n-σ-D为一个大小为M×N的矩阵，N表示所述单分量立方相位信号的时域长度，M表示所述单分量立方相位信号的频域长度；参数空间n-σ-D第l行第k列的取值为 nk表示离散时间变量n的第k个取值，k＝1,2,...,N，N为单分量立方相位信号的长度；σl表示离散频率变量σ的第l个取值，l＝1,2,...,M。
[0055] 从步骤3中得出，参数α和β的估计问题转换为直线检测问题；参照图2a，为本发明实施例中直线检测的第一原理示意图，参照图2b，为本发明实施例中直线检测的第一原理示意图；图2a中，横轴表示参数n，纵轴表示参数σ；图2b中，横轴表示参数β，纵轴表示参数α。只需求取直线σ-α-βn＝0的截距和斜率便可获得参数α和β的估计值。根据直线检测原理图知，在以参数α为截距、β为斜率的直线σ＝α+βn所在的平面n-σ上，直线σ＝α+βn上的每一点在平面β-α上都对应于一条直线，直线σ＝α+βn上的每一点在平面β-α上对应的这些直线相交于同一点。
[0056] 本发明实施例中，经验性地设截距α共有M个离散化的取值，分别表示为α1至αM；经验性地设斜率β共有N个离散化的取值，分别表示为β1至βN。根据上述分析，建立参数空间β-α-Ψ，参数空间β-α-Ψ中两个轴的量化尺度与原参数空间n-σ-D的量化尺度相同，即同样是一个大小为M×N的矩阵；令i＝1,2,...,N，令j＝1,2,...,M，参数空间β-α-Ψ第j行第i列的取值 应为：满足方程σ＝αj+βin的所有二维坐标(n,σ)对应的时频分布的模值之和，αj表示截距α的第j个取值，βi表示斜率β的第i个取值。可以看出，参数空间β-α-Ψ第i列可以由参数空间n-σ-D每一列进行循环移位累加而成。下面对循环移位过程进行具体说明。
[0057] 令i＝1,2,...,N，当βink小于0时，将参数空间n-σ-D第k列数据向下循环移位-βink位，得出参数空间n-σ-D第k列数据的第i循环移位数据；当βink大于等于0时，将参数空间n-σ-D第k列数据向上循环移位βink位，得出参数空间n-σ-D第k列数据的第i循环移位数据。
[0058] 根据上述分析可知，参数空间n-σ-D第k列数据是大小为M×1的列向量，参数空间n-σ-D第k列数据可以表示为：
[0059]
[0060] 其中，上标T表示矩阵或向量的转置。
[0061] 在将参数空间n-σ-D第k列数据向下循环移位或向上循环移位之前，得出|βink|MOD M的值，MOD表示求余运算；如果|βink|MOD M＝0，则参数空间n-σ-D第k列数据向上循环移位0位，即参数空间n-σ-D第k列数据不进行循环移位，此时参数空间n-σ-D第k列数据的第i循环移位数据为参数空间n-σ-D第k列数据。如果|βink|MOD M≠0且βink小于0，将参数空间n-σ-D第k列数据向下循环移位|βink|MOD M位，得出参数空间n-σ-D第k列数据的第i循环移位数据；如果|βink|MOD M≠0且βink大于0，将参数空间n-σ-D第k列数据向上循环移位|βink|MOD M位。
[0062] 当βink小于0时，参数空间n-σ-D第k列数据的第i循环移位数据是大小为M×1的列向量，参数空间n-σ-D第k列数据的第i循环移位数据的第n1行元素为参数空间n-σ-D第k列数据的第|βink|MOD M+n1行元素，n1取1至M-(|βink|MOD M)；参数空间n-σ-D第k列数据的第i循环移位数据的第n2行元素为参数空间n-σ-D第k列数据的参数空间n-σ-D第k列数据的第n2-M+(|βink|MOD M)行元素，n2取M-(|βink|MOD M)+1至M。
[0063] 当βink大于0时，参数空间n-σ-D第k列数据的第i循环移位数据是大小为M×1的列向量，参数空间n-σ-D第k列数据的第i循环移位数据的第n3行元素为参数空间n-σ-D第k列数据的第M-(|βink|MOD M)+n3行元素，n3取1至(|βink|MOD M)；参数空间n-σ-D第k列数据的第i循环移位数据的第n4行元素为参数空间n-σ-D第k列数据的参数空间n-σ-D第k列数据的第n4-(|βink|MOD M)行元素，n4取(|βink|MOD M)+1至M。
[0064] 特别地，当βink等于0时，参数空间n-σ-D第k列数据向上循环移位0位，即参数空间n-σ-D第k列数据不进行循环移位，此时参数空间n-σ-D第k列数据的第i循环移位数据为参数空间n-σ-D第k列数据。
[0065] 下面以几个具体的例子对本发明实施例的循环移位过程进行说明。
[0066] 设M＝9，则参数空间n-σ-D第k列数据可以表示为：
[0067]
[0068] 如果βink为0或βink为9的整数倍，参数空间n-σ-D第k列数据的第i循环移位数据为参数空间n-σ-D第k列数据；如果βink＝3，则|βink|MOD M＝3，此时，参数空间n-σ-D第k列数据的第i循环移位数据为：
[0069]
[0070] 如果βink＝-30，则|βink|MOD M＝3，此时，参数空间n-σ-D第k列数据的第i循环移位数据为：
[0071]
[0072] 本发明实施例中，在得出参数空间n-σ-D第k列数据的第i循环移位数据之后，将参数空间n-σ-D第1列数据的第i循环移位数据至参数空间n-σ-D第N列数据的第i循环移位数据进行累加，得出参数空间β-α-Ψ的第i列 是大小为M×1的列向量；则参数空间β-α-Ψ为：
[0073]
[0074] 步骤5，利用峰值检测检测出参数空间β-α-Ψ元素最大值对应的行序号j(max)和参数空间β-α-Ψ元素最大值对应的列序号i(max)；得出上述实际生活中获取的信号的调频率的估计值 以及上述实际生活中获取的信号的二次调频率的估计值 为截距α的第j(max)个取值αj(max)， 为斜率β的第i(max)个取值βi(max)。
[0075] 在得出上述实际生活中获取的信号的调频率的估计值 以及述实际生活中获取的信号的二次调频率的估计值 之后，利用得到估计值构造补偿函数
用该补偿函数来补偿原始信号的高阶相位项，得到新信
号 通过傅里叶变换完成补偿后信号的
能量的积累，对能够积累后的信号进行峰值检测，得出信号幅度A的估计值以及所述单分量立方相位信号的中心频率f的估计值。
[0076] 下面结合仿真实验对本发明效果做进一步的说明和分析
[0077] 1)仿真条件：
[0078] 本发明的仿真的立方相位参数为信号幅度A＝1，中心频率f＝20Hz，调频率α＝40Hz/s，二次调频率β＝60Hz/s2。信号长度N＝512，采样频率Fs＝256Hz。
[0079] 2.仿真内容：
[0080] 按照上述仿真条件，在
软件MATLAB8.0中进行仿真实验，具体如下：
[0081] 求离散单分量立方相位信号的瞬时自相关函数，并沿着时延轴做NUFFT，得到该单分量立方相位信号的时频分布，参照图3，为仿真实验中得出的单分量立方相位信号的时频分布的示意图；图3中，
水平面的两个轴分别表示时间和频率，竖直轴表示时频分布的归一化模值(能量)；参照图4，为图3中所示的时频分布的等高线图。图4中，横轴表示时间，纵轴表示频率。从图4中可以清晰地出，时间-频率是呈一条直线分布的，即能量集中于该直线，且该直线的截距和斜率分别为立方相位信号的调频率α和二次调频率β。
[0082] 对该单分量立方相位信号的时频分布进行循环移位操作，参照图5，为仿真实验中得出新的参数空间的示意图。在图5中，变换
坐标系，将整个n-σ-D数据空间的所有数据，按照顺序移位原则进行循环移位并累加，得到新的参数空间β-α-Ψ。经过上述算法计算后得到二次调频率-调频率分布，参照图6，为仿真实验中得出的二次调频率-调频率分布示意图。图6中，水平的两个轴分别表示二次调频率和调频率，竖直轴表示归一化幅度(能量)，参照图7，为图6的等高线示意图。图7中，横轴表示二次调频率，纵轴表示调频率。
[0083] 3)仿真结果分析：
[0084] 从图2可以看出，通过对把直线上点的坐标变换到过点的直线的系数域，便可以将参数估计问题转化为变换域直线检测问题。
[0085] 从图3、图4可以看出，通过对信号自相关函数进行NUFFT操作后，立方相位信号的能量主要集中于一条直线，且该直线的截距和斜率分别为立方相位信号的调频率α和二次调频率β。因此，可以通过直线检测方法实现立方相位信号的相位参数估计。
[0086] 从图5可以看出，先建立新的参数空间后，再利用循环移位累加方法，可以十分便捷地完成两个参数空间的转换。
[0087] 从图6、7可以看出，使用本发明所提出的方法后，仿真实验所得到的结果与之前仿真条件中所设定的参数完全一致。由此证明，本发明的方法确实有效。
[0088] 综上，本发明降低了立方相位信号参数估计的复杂度和计算量，避免了复杂的二维谱峰搜索过程，提高了估计精度和处理效率。
[0089] 显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些
修改和变型属于本发明
权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。
