贝塞尔函数与虚拟阵列相结合的SFM信号的参数估计方法

摘要

本发明请求保护贝塞尔函数和虚拟阵列相结合的SFM信号的参数估计方法，属于信号处理技术领域。本方法充分考虑到SFM信号的复杂性及非平稳特性，采用了一种特殊的处理方法，把其变换成那些常用的窄带信号形式，根据分解后的窄带信号的特点，对其做相应的处理，然后利用估计信号方向的方法来估计该信号的频率参数；其次利用超定方程组最小二乘解的定理解出其调制系数。该方法不仅能够很好的处理SFM信号，而且能够准确的估计出信号的参数，同时可以推广到多分量的信号。

说明书

技术领域

本发明涉及信号与信息处理领域，具体为一种复SFM（正弦调频）信号 的参数估计方法。

背景技术

复正弦调频信号是一种典型的非线性调频信号，是一种具有低截获率的 时变信号，具有抑制泄露和近区干扰等特性，因此广泛应用在雷达、声纳、 通信、生物工程以及地震勘测等领域。

目前，人们已经提出了多种非线性调频信号形式，主要包括多项式相位 信号和复SFM信号。虽然非线性调频信号的参数估计有一系列的方法，但是 大多数都是针对多项式相位信号的，而对复SFM信号的研究较少。文献 (Braham Barket.Instantaneous frequency estimation of nonlinear  frequency-modulated signals in the presence of multiplicative additive noise[J] IEEE Trans.on SP,2001,49(10):2214-2222.)充分讨论了在不同信噪比下的非 线性调频信号在多项式Wigner-Ville分布（PWVD）平面上瞬时频率的反映特 征，但是也是把非线性调频信号统一建模成多项式相位信号进行讨论的。故 很多方法只能处理多项式相位信号，而不能处理复SFM信号。文献（熊刚， 赵惠昌，王李军.伪码-载波调制侦察信号识别的谱相关方法（Ⅱ）--伪码-载 波调频信号的调制识别和参数估计[J].电子与信息学报.2005,27(7): 1087-1092.）提出了利用谱相关方法对伪码SFM信号进行调制识别和参数估 计，但是该算法运算量庞大，并且只适用于窄带SFM信号。文献（熊辉，吕 远，增德国，唐斌.利用卡森准则的正弦调频信号参数估计方法[J].电子测量 与仪器学报,2010,4(24):353-358.）利用信号频谱对称特性，给出了信号载频 估计的实现方法，基于瞬时频率的频谱分析，推导了调制频率的估计表达式， 但是对于多个复SFM信号，该参数估计算法会失效，并且在低信噪比情况下 参数估计性能欠佳。

然而在所有的参数估计算法中，很多常规的算法在理想条件下具有良好 的性能，但是当信号源相干或者是宽带的信号时性能就会变得很差，甚至无 法进行相应的参数估计。综上所述，对复SFM信号的参数估计研究较少，现 有算法不完善，或受带宽限制，或只能估计部分参数。

发明内容

本发明所要解决的技术问题是，提出一种贝塞尔函数与虚拟阵列相结合 的SFM信号的参数估计方法，该方法克服了传统方法只能处理窄带信号的困 难，解决了SFM信号参数估计性能差的问题。在信号处理过程中，可以有效 地把宽带的SFM信号分解成一组具有谐波幅度、谐波频率的窄带信号模型的 形式，然后结合虚拟阵列变换法进行解相干，从而利用常规的子空间方法进 行参数估计。该方法不仅计算量小，而且估计性能明显优于未经解相干时的 参数估计。

本发明解决上述技术问题的技术方案是，在贝塞尔函数及雅可比变换处 理复正弦调频信号的基础上，采用虚拟阵列变换法对分解后的窄带正弦信号 进行解相干，从而利用信号子空间（MUSIC（多重信号分类）算法等）的方 法得到信号谐波频率估计，根据这些谐波频率的对称性，从而得到载波频率 和调制频率估计。SFM信号参数估计方法的步骤具体包括：

给定信号s(t)，在t=nΔt时刻对给定信号s(t)离散化得到离散采样序列为 {s(n)}，其中t是时间变量，Δt为采样间隔。对离散化后的信号按照雅可比展 开恒等式：进行展开。其中，J_v(k)是第一 类v阶贝塞尔函数，j为虚数单位。展开后可得无穷多个具有谐波幅度、谐波 频率的正弦信号之和的形式；将接收信号x(n)=s(n)+W(n)中s(n)部分分量 忽略不计，其中W(n)代表高斯白噪声；对忽略s(n)分量后的接收信号进行分 解，把其分解成窄带信号模型为：X=BC+W。其中，B为方向向量，C为 信号矢量，W为高斯白噪声。利用虚拟阵列变换法对分解后的窄带信号进行 解相干；具体可为：对真实阵列，采样得信号1，调用公式： R_l=B_lG^l-1R_ls(G^l-1)^H(B_l)^H+W求得真实阵列的协方差R₁，其中G^l表示Q×Q 对角阵G的l次方：Q表示B中方向向量个 数，d表示虚拟阵列平移的间距，f_i表示对应的第i个分量的谐波频率，R_ls是 信号的协方差矩阵。然后利用内插技术得第l个虚拟阵列的协方差R_l，调用公 式：其中，表示信号协方差矩阵的平均，求 得平均协方差矩阵利用常规方法（如MUSIC（多重信号分类）算法等） 估计出各个谐波频率。根据各个谐波频率的对称性质可以得到各个谐波频率 之和正好是载波频率的整数倍，计算任意两个相邻谐波频率之差获得调制 频率根据第一类贝塞尔函数的最高阶数的递归关系： J_v+1(m_f)=2vJ_v(m_f)/m_f-J_v-1(m_f)，（式中J_k(·)为第一类k阶Bessel函数，k为 Bessel函数的阶数），组成超定方程组。求得方程组最小二乘解，从而估计出 其调制系数

该方法充分利用了虚拟阵列变换法和信号子空间方法，将变换后的窄带 信号有效地的进行了解相干，降低了参数估计的误差，极大的减少了算法的 计算量，提高了参数估计的精度。相比已有的参数估计方法从计算量、估计 性能方面得到显著的提高。

附图说明

图1本发明贝塞尔函数和虚拟阵列结合方法的处理框图；

图2本发明SFM信号参数估计流程图；

图3虚拟阵列变换法原理图；

图4m_f=1时的谐波幅度；

图5m_f=2时的谐波幅度；

图6本发明变换后的窄带信号实部；

图7变换后各个谐波分量的谱线图；

图8本发明解相干后MUSIC算法的谐波频率估计；

图9本发明未经解相干的MUSIC算法谐波频率估计；

图10本发明解相干后谐波频率估计的误差分析。

具体实施方式

下面结合附图和实例，对本发明的实施作进一步详细说明，但本发明的 实施方式并不仅限于此。

如图1所示为本发明贝塞尔函数和虚拟阵列结合的参数估计方法示意图， 首先对接收信号进行离散化，根据雅可比变换及贝塞尔函数的性质对其进行 窄带化处理，然后把变换后的信号写成窄带信号模型的形式，根据该信号的 特点采用虚拟阵列变换法进行解相干，然后采用信号子空间的方法进行谐波 频率的估计，以及根据贝塞尔函数的递归性质，组成超定方程组，从而解得 其调制系数。

如前所述，本发明是在贝塞尔函数及雅可比变换处理复SFM信号的基础 上对SFM信号进行参数估计的。复SFM信号表示为：

s(t)=Aexp{j[2πf_ct+m_fsin(2πf_mt)]},0≤t≤T    (1)

其中，A为信号的幅度，f_c为信号载波频率，f_m为信号调制频率，m_f为调频 系数，T为信号周期。

在t=n△t时刻对给定信号s(t)离散化得到离散采样序列为{s(n)}：

s(n)=Aexp{j{2πf_c(nΔt)+m_fsin[2πf_m(nΔt)]}},0≤n≤N-1    (2)

其中，N为采样长度，Δt为采样间隔。

接收信号可表示为：

x(n)=s(n)+W(n)

=Aexp{j{2πf_c(nΔt)+m_fsin[2πf_m(nΔt)]}}+W(nΔt),0≤n≤N-1    (3)

其中，W为零均值的高斯白噪声。

对离散化后的信号根据雅可比展开恒等式：

$\exp \left(\mathrm{jk}\sin \beta \right)=\underset{v=-\infty }{\overset{\infty }{\Sigma }}{J}_v\left(k\right)\exp \left(\mathrm{jv\beta }\right)---\left(4\right)$

其中，J_v(k)是第一类v阶贝塞尔函数。把式(4)代入式(3)可得无穷多个具有谐 波幅度谐波频率的正弦信号之和的形式，一般情况下取Δt=1，即：

$x\left(n\right)=s\left(n\right)+W\left(n\right)$

$=\mathrm{Aexp}\left\{j\right[2\pi f_{c}n+m_f\sin \left(2\pi f_{m}n\right)\left]\right\}+W\left(n\right)$

$=\underset{v=-\infty }{\overset{\infty }{\Sigma }}{\mathrm{AJ}}_v\left(m_f\right)\exp \left[j2\pi (f_c+{\mathrm{vf}}_m)n\right]+W\left(n\right),0\le n\le N-1---\left(5\right)$

由式(5)可以看出，非线性的复正弦调频信号被变成了线性的正弦信号形式。

其中，由贝塞尔函数的对称性质： J_v(m_f)=(-1)^vJ_-v(m_f)，并且当|v|>|m_f|时，J_v(m_f)≈0，s(n)部分分量小到 可以被忽略不计。故式(5)可以写成：

$x\left(n\right)\approx \underset{v=-V}{\overset{V}{\Sigma }}{\mathrm{AJ}}_v\left(m_f\right)\exp \left[j2\pi (f_c+{\mathrm{vf}}_m)n\right]+W\left(n\right),0\le n\le N-1---\left(6\right)$

其中，V是信号分量不可以被忽略的第一类贝塞尔函数的最高阶数。当m_f>1 时，V≈m_f+1，然而当m_f∈[0.14,1]时，V=1或2；当m_f∈[00.14)时，V=0。 经过分解，式(6)有下面的表达形式：

$x\left(n\right)\approx \underset{v=-V}{\overset{V}{\Sigma }}{\mathrm{AJ}}_v\left(m_f\right)\exp \left[j2\pi (f_c+{\mathrm{vf}}_m)n\right]+W\left(n\right),0\le n\le N-1$

$={\mathrm{AJ}}_{-V}\left(m_f\right)\exp \left(j2\pi (f_c-{\mathrm{Vf}}_m)n\right)+···+{\mathrm{AJ}}_0\left(m_f\right)\exp \left(j2\pi (f_c+0f_m)n\right)$

$+···+{\mathrm{AJ}}_v\left(m_f\right)\exp \left(j2\pi (f_c+{\mathrm{Vf}}_m)n\right)+W\left(n\right)$

$=\mathrm{BC}+W$

$=[\exp \left(j2\pi (f_c-{\mathrm{vf}}_m)n\right)···\exp \left(j2\pi (f_c+{0f}_m)n\right)···\exp \left(j2\pi (f_c+{\mathrm{Vf}}_m)n\right)]·$

${[{\mathrm{AJ}}_{-V}\left(m_f\right)···{\mathrm{AJ}}_0\left(m_f\right)···{\mathrm{AJ}}_V\left(m_f\right)]}^T+W\left(n\right)$

(7)

其中，B=[exp(j2π(f_c-Vf_m)n)…exp(j2π(f_c+0f_m)n)…exp(j2π(f_c+Vf_m)n)]， C=[AJ_-V(m_f)…AJ₀(m_f)…AJ_V(m_f)]^T。

其中，f₁=f_c-Vf_m,…,f_V+1=f_c+0f_m,…,f_2V+1=f_c+Vf_m为各个谐波频率。由式(7) 可以看出接收的信号可以表示成窄带信号模型的形式，即X=BC+W，在此B 相当于方向向量，C相当于信号矢量，W为噪声。由 知，C中的几个分量都是常数，故相当于每个信 号只是相差常数倍，即如果把它们看做是信号矢量的话，它们是相干的。

贝塞尔函数方法能够把难处理的宽带复正弦调频信号分解成一组具有谐 波幅度、谐波频率的窄带信号形式。

由(7)式可以看出，复正弦调频信号经过贝塞尔函数及雅可比展开式可以 被分解成窄带信号模型的形式。然而其信号矢量却是相干的，因为相干的信 号可以合并成一个信号，导致信号子空间的维数小于信号源的数量。信号的 自相关矩阵成为非满秩矩阵，因此也就不能进行正确的估计。为解决这个难 题，本发明利用虚拟阵列变换的方法，提出了贝塞尔函数与虚拟阵列变换结 合的参数估计方法。

虚拟阵列变换法就是每次将M个阵元的均匀线阵等间距的向右平移一个 等距d，形成一系列的虚拟阵列，如图3所示。对于其中的每一个接收序列， 第l个子阵的接收的数据矢量表示为X_l(t)，且有：

X_l(t)=B_lG^l-1C+W_l          (8)

整个阵列等间距地向右平移一个d，则平移l-1次后，得第l(1≤l≤M′)个阵 列协方差矩阵为：

R_l=B_lG^l-1R_ls(G^l-1)^H(B_l)^H+σ²I      (9)

其中，$B_l=[b_l\left({\theta }_1\right),b_l\left({\theta }_2\right),···b_l\left({\theta }_Q\right)],b\left({\theta }_i\right)={[1,e^{-j2\pi d\sin {\theta }_i/\lambda },···,e^{-j2\pi (l-1)d\sin {\theta }_i/\lambda }]}^T,$另外， β_i=(2πdsinθ_i)/λ,i=1,2,…,Q，Q=2V+1为信号源数，d为均匀线阵的间距，c 为信号传播速度。G^l表示Q×Q对角阵G的l次方。则：

$G=\mathrm{diag}[e^{j2\pi {\mathrm{df}}_1},e^{j2\mathrm{\pi d}{f}_2},···,e^{j2\pi {\mathrm{df}}_Q}]---\left(10\right)$

其中，f₁,f₂,…f_Q表示各个谐波频率。

则整个线阵的协方差矩阵定义为各子阵协方差的平均，即：

$\bar{R}=\frac{1}{M^{\prime }}\underset{l=1}{\overset{M^{\prime }}{\Sigma }}{R}_l=B_l{\bar{R}}_{\mathrm{ls}}{B}_l^H+{\sigma }^{2}I---\left(11\right)$

${\bar{R}}_{\mathrm{ls}}=\frac{1}{M^{\prime }}\underset{l=1}{\overset{M^{\prime }}{\Sigma }}{G}^{l-1}{R}_{\mathrm{ls}}{G}^{-(l-1)}---\left(12\right)$

如果是满秩的，就可以利用它进行相干信号的参数估计。由以上的分 析可以Q个信号是完全相干的，显然信号的协方差矩阵R_ls的秩为1，即：

R_ls=ββ^H              (13)

其中，是一个行向量，β₁β₂…β_Q为其各个分量。于是 有：

${\bar{R}}_{\mathrm{ls}}=\frac{1}{M^{\prime }}\underset{l=1}{\overset{M^{\prime }}{\Sigma }}{G}^{(l-1)}\beta {\beta }^H{G}^{-(l-1)}$

$=\frac{1}{M^{\prime }}\underset{l=1}{\overset{M^{\prime }}{\Sigma }}\left[G^{(l-1)}\beta \right]{\left[G^{(l-1)}\beta \right]}^H$

$=\frac{1}{M^{\prime }}\mathrm{DD}---\left(14\right)$

其中，D=[β,Gβ,…,G^(l-1)β]。显然有的秩和D的秩相等，矩阵D可表示为：

$D=\left(\begin{array}{cccc}{\beta }_1{e}^{j2\pi {\mathrm{df}}_1}& {\beta }_1{e}^{j2*2\pi {\mathrm{df}}_1}& ···& {\beta }_1{e}^{j2\pi M^{\prime }{\mathrm{df}}_1}\\ {\beta }_2{e}^{j2\pi {\mathrm{df}}_2}& {\beta }_2{e}^{j2*2\pi {\mathrm{df}}_2}& ···& {\beta }_2{e}^{j2\pi M^{\prime }{\mathrm{df}}_2}\\ & ·& ·& ·\\ & ·& ·& ·\\ & ·& ·& ·\\ {\beta }_Q{e}^{j2\pi {\mathrm{df}}_Q}& {\beta }_Q{e}^{j2*2\pi {\mathrm{df}}_Q}& ···& {\beta }_Q{e}^{j2\pi M^{\prime }{\mathrm{df}}_Q}\end{array}\right)---\left(15\right)$

从式(15)可以看出，只要所有的信号不是来自同一个方向，矩阵D满秩。因而 矩阵和矩阵都为满秩的，即：

$\mathrm{Rank}\left[{\bar{R}}_{\mathrm{ls}}\right]=\mathrm{Rank}\left[D\right]=\min (Q,M^{\prime })---\left(16\right)$

只要满足M′≥Q，就有于是恢复了满秩，由前面的分 析可知，SFM信号被分解成了窄带信号模型的形式，解相干后从而可以使矩 阵对相干信号用常规的子空间分解方法（MUSIC算法）来进行估计。为了 使阵列的自由度得到充分利用，一般取M′=M，即子阵列的个数与阵列的阵 元数相等。在此由于没有具体的接收阵列，故让子阵列的个数、阵列的阵元 数、信号的采样长度都相同。

根据贝塞尔函数的对称性，当各个谐波频率相加时，其和正好是载波频 率的整数倍，任意两个相邻谐波频率之差正好是调制频率，从而估计出载波 频率和调制频率。以上是信号的载波频率、调制频率的估计。

为了解决传统方法的局限性，受参数限制或者只能估计部分参数的问题， 本发明还可采用调制系数的估计。

由第一类贝塞尔函数，可以得到一下递归关系：

J_v+1(m_f)=2vJ_v(m_f)/m_f-J_v-1(m_f)        (17)

由式(17)的递归关系，假设贝塞尔函数分解的最高阶数为V，则贝塞尔函数的 各个阶数分别为-V-V+1…V。则有下列式成立：

$\left(\begin{array}{cc}v=V-1& J_V\left(m_f\right)=2(V-1)J_{V-1}\left(m_f\right)/m_f-J_{V-2}\left(m_f\right)\\ v=V-2& J_{V-2}\left(m_f\right)=2(V-2)J_{V-2}\left(m_f\right)/m_f-J_{V-3}\left(m_f\right)\\ ···& ···\\ v=-V+1& J_{-V+2}\left(m_f\right)=2(-V+1)J_{-V+1}\left(m_f\right)/m_f-J_{-V}\left(m_f\right)\end{array}\right)---\left(18\right)$

而式(18)又可以分解成如下形式：

$J=\left(\begin{array}{c}{J}_V\left(m_f\right)\\ J_{V-2}\left(m_f\right)\\ ·\\ ·\\ ·\\ J_{-V+2}\left(m_f\right)\end{array}\right)=\left(\begin{array}{cc}2(V-1)J_{V-1}\left(m_f\right)& -J_{V-2}\left(m_f\right)\\ 2(V-2)J_{V-2}\left(m_f\right)& -J_{V-3}\left(m_f\right)\\ ·& ·\\ ·& ·\\ ·& ·\\ ·& ·\\ ·& ·\\ ·& ·\\ 2(-V+1)J_{-V+1}\left(m_f\right)& -J_{-V}\left(m_f\right)\end{array}\right)·\left(\begin{array}{c}1/m_f\\ 1\end{array}\right)=\mathrm{Zb}---\left(19\right)$

根据超定方程组的定义可知，式(19)是超定方程组，所谓超定方程组是指 方程个数大于未知量个数的方程组。对超定方程组：

当Z^TZ可逆时，超定方程组(19)存在最小二乘解，且即为方程组 Z^TZb=Z^TJ的解，即b=(Z^TZ)^-1Z^TJ。于是有：$\left(\begin{array}{c}1/m_f\\ 1\end{array}\right)={\left(Z^{T}Z\right)}^{-1}{Z}^{T}J.$

如图2所示为对SFM信号参数估计流程示意图，具体步骤如下，假设有 两个相同时频分布的宽带SFM信号：

s₁(t)=A₁exp{j[2π*f_ct+m_f*sin(2π*f_mt)]}

s₂(t)=A₂exp{j[2π*f_ct+m_f*sin(2π*f_mt)]}

信号参数设置为：m_f=1，f_c=0.35，f_m=0.03，信号幅度为A₁=1，A₂=2A₁， SNR＝20dB，快拍数N=256。

步骤1：在t=nΔt时刻对给定信号s₁(t)、s₂(t)离散化得到离散采样序列 为{s₁(n),s₂(n)}。

步骤2：利用雅可比展开式：并根据 贝塞尔函数变换后得到窄带信号x(n)可以表示为：

$x\left(n\right)\approx \underset{v=-2}{\overset{2}{\Sigma }}{\mathrm{AJ}}_v\left(m_f\right)\exp \left[j2\pi (f_c+{\mathrm{vf}}_m)n\right]+W\left(n\right),0\le n\le N-1---\left(20\right)$

经过分解，式(20)可以写成下面的表达形式：

$x\left(n\right)\approx \underset{v=-2}{\overset{2}{\Sigma }}{\mathrm{AJ}}_v\left(m_f\right)\exp \left[j2\pi (f_c+{\mathrm{vf}}_m)n\right]+W\left(n\right),0\le n\le N-1$

$=A{J}_{-2}\left(m_f\right)\exp \left[j2\pi (f_c-2f_m)n\right]+A{J}_{-1}\left(m_f\right)\exp \left[j2\pi (f_c-1f_m)n\right]+A{J}_0\left(m_f\right)\exp \left[j2\pi (f_c+0f_m)n\right]$

$+{\mathrm{AJ}}_1\left(m_f\right)\exp \left[j2\pi (f_c+1f_m)n\right]+{\mathrm{AJ}}_2\left(m_f\right)\exp \left[j2\pi (f_c+2f_m)n\right]+W\left(n\right)$

$=\mathrm{BC}+W---\left(21\right)$

步骤3：式(21)虽然写成了窄带信号模型的形式，但是由于由各个贝塞尔 函数组成的信号是相干的，故不能直接利用常规的窄带信号算法对其进行估 计，在此利用虚拟阵列变化法对其解相干，从而利用基于信号子空间的MUSIC 算法等对其进行谐波频率的估计。

步骤4：根据各个谐波频率的特点从而估计出载波频率和调制频率。

步骤5：由式(17)、式(18)可知，当m_f=1时，即第一类贝塞尔函数的最高 阶数为2时，该递归关系可以表示为：

$\left(\begin{array}{cc}v=1,& J_2\left(m_f\right)={2J}_1\left(m_f\right)/m_f-J_0\left(m_f\right)\\ v=0,& J_1\left(m_f\right)=2·0·J_0\left(m_f\right)/m_f-J_{-1}\left(m_f\right)\\ v=-1,& J_0\left(m_f\right)=2·(-1)·J_{-1}\left(m_f\right)/m_f-J_{-2}\left(m_f\right)\end{array}\right)---\left(22\right)$

于是由式(22)可以得到如下关系：

$J=\left(\begin{array}{cc}{2J}_1\left(m_f\right)& -J_0\left(m_f\right)\\ 2·0·J_0\left(m_f\right)& -J_{-1}\left(m_f\right)\\ 2·(-1)·J_{-1}\left(m_f\right)& -J_{-2}\left(m_f\right)\end{array}\right)·\left(\begin{array}{c}\frac{1}{m_f}\\ 1\end{array}\right)=J_Z\left(\begin{array}{c}\frac{1}{m_f}\\ 1\end{array}\right)---\left(23\right)$

故由$\left(\begin{array}{c}\frac{1}{m_f}\\ 1\end{array}\right)={\left(J_Z^T{J}_Z\right)}^{-1}$J从而估计出调制系数

为了更客观的说明该方法的有效性。图4、图5是本发明不同调制频率的 谐波幅度图，该组图证明了第一类贝塞尔函数的对称性质以及v≈m_f+1，当 |v|>|m_f|时，J_v(m_f)≈0，即s(n)小到可以被忽略不计。图4取m_f=1，当V≥3 时，可以明显看出J₃(1)≈0；图5取m_f=2，当V≥4时，可以明显看出J₄(2)≈0。

由式(7)、(21)可知，被贝塞尔函数分解后的信号变成了窄带信号模型的形 式，图6给出了变换后的窄带信号分别在没有噪声和SNR=20dB时，信号的实 部幅值特性图。从图中可以看出经过变换后，当没有噪声时，信号变换为窄 带的信号；当SNR=20dB时，接收的信号变换为窄带信号和噪声。从而把宽带 的信号变成了窄带的信号，有利于用比较成熟的方法进行参数估计。图7给 出了变换后各个谐波分量的谱线图。即与图4中的各个非零部分的谐波频率 相对应。

由以上分析可知，SFM信号经过贝塞尔函数变换及虚拟阵列解相干后， 可以利用常规的基于子空间的MUSIC算法等对其进行谐波频率的估计，本发 明采用MUSIC算法估计其谐波频率。为了充分说明贝塞尔函数和虚拟阵列相 结合的有效性，本发明做了一组对比试验，即贝塞尔函数不和虚拟阵列变换 法结合直接进行谐波频率的估计。其仿真结果见图8、图9所示。在前面的信 号参数已知的情况下，各个谐波频率的理论值应该为：f₁=(f_c-2f_m)=0.29， f₂=(f_c-1f_m)=0.32，f₃=(f_c+0f_m)=0.35，f₄=(f_c+1f_m)=0.38， f₅=(f_c+2f_m)=0.41。

由图8、图9可以看出，该方法能够精确的估计出各个谐波频率。通过对 比分析，解相干后的MUSIC算法估计出来的频率谱峰更加尖锐，更加突出； 而未经解相干的频率估计，明显没有解相干后的频率估计要好，特别是其边 上的两个频率估计几乎看不到谱峰，从而证明了贝塞尔函数和虚拟阵列变换 法相结合的优越性。然而由图8估计出各个谐波频率：f₁=(f_c-2f_m)， f₂=(f_c-1f_m)，f₃=(f_c+0f_m)，f₄=(f_c+1f_m)，f₅=(f_c+2f_m)。然后再根据 这些谐波频率的对称性可以看出，所有谐波频率之和是载波频率的整数倍， 即任意两个相邻谐波频率之差是调制频率，即 或者其中载波频率的真实值为f_c=0.35，仿真估计得 到的是而调制频率真实值为f_m=0.03，仿真估计得到的是 因此在误差允许范围内能够精确的估计出载波频率和调制频率。

然后利用贝塞尔函数的递归性质，组成超定方程组，根据超定方程组估 计出其调频系数。为了充分的证明该方法的有效性，进一步实验，让调制系 数分别为m_f=1和m_f=2。

根据超定方程组的定理可知，估计出来的应该是调制频率的倒数，因此 估计出的结果为：由仿真结果可知，估计出来的调频 系数和真实值基本上是完全相同，因此调频系数得到了很好的估计。

为了说明谐波频率的估计性能，本发明做了一组实验是各个谐波频率的 误差分析，其仿真结果见图10。从图10中可以看出，各个谐波频率的标准差 随着信噪比的增加在逐渐的降低，也即说明了其估计性能越来越好。

本发明采用的是具有相同时频分布的2个分量的SFM信号，该方法很容 易被推广到多个分量的SFM信号的情形。并且该方法不仅能够精确的估计出 各个参数，而是不受部分参数限制等。