ISSN 
下载:
-
多音信号的频率估计是许多领域的基本问题，包括雷达[1]和无线通信[2]。如在正交频分复用(orthogonal frequency division multiplexing, OFDM)通信系统的载波频偏估计[2]和多音干扰消除[3]中会遇到多音信号的频率估计。现有的用于多音信号频率估计算法主要分为两类：时域估计算法和频域估计算法。
时域估计算法通常会需要矩阵求逆、特征值分解、奇异值分解[4]，因此其计算复杂度为
$O( {L{N^3}} )$ ，其中$L$ 为单音分量的数量，$N$ 为采样点数。频域估计算法主要基于离散傅里叶变换(discrete fourier transform, DFT)和一些插值形式，因此他们的计算复杂度一般为
$O\left( {LN{{\log }_2}N} \right)$ ，这远小于时域估计算法，因此频域估计算法得到了更多的关注和应用。针对单音信号，已有很多插值器被提出[5]，但在处理多音信号时，频谱泄露造成的分量间互干扰会降低这些插值器的估计精度。为了降低分量间互干扰，许多基于加窗的插值器被提出。文献[6]提出基于赖夫−文森特类I窗的插值器。文献[7]提出汉宁窗，文献[8]提出余弦窗。文献[9]提出了一种适用于经典窗函数的补零插值器，然而由于其中的三次方程，该插值器的估计性能较差。文献[10]提出了一种高阶多项式插值算法来补偿特定于窗函数的估计器偏差，支持任意窗函数。基于文献[10]，文献[11]提出了一个新的插值器，该插值器仅根据窗函数调整偏差补偿因子，因此其复杂度更低。虽然这些非矩形窗可以减少分量间的互干扰，但也会造成一定的功率损失。文献[12]提出了牛顿化正交匹配追踪(newtonalized orthogonal matching pursuit, NOMP)算法，该算法通过牛顿迭代法和反馈提高频率估计的精度，然而其计算复杂度较高。文献[13]提出了一种迭代估计算法，该算法采用了A&M插值器[14]。文献[15]提出了一种两阶段估计算法，该算法采用了抛物线插值器[16]来克服分量间互干扰，然而当分量间频率间隔较小时，其估计性能仍然较差。为了降低分量间的相互干扰，提高频率估计性能，本文提出了一种基于两阶段加窗插值(two-stage windowed interpolation, TSWI)的多音信号频率估计算法，该算法包含一个新的支持任意窗函数的插值器。在第一个阶段，通过选用可以降低频谱泄露的窗函数获得初始估计值；在第二个阶段，由于分量间的相互干扰可以通过初始估计值消除，通过选用矩形窗函数弥补第一阶段非矩形窗带来的信噪比(signal-to-noise ratio, SNR)损失并得到精估计值。数值仿真结果表明，本文算法具有比现有算法更优的估计性能。
-
含噪声的多音信号可以表示为：
$$r\left[ n \right] = \sum\limits_{l = 1}^L {{s_l}\left[ n \right]} + w\left[ n \right] = \sum\limits_{l = 1}^L {{A_l}{{\rm{e}}^{{\rm{j}}2{\text π} \frac{{{f_l}}}{{{f_s}}}n}}} + w\left[ n \right]$$ (1) 式中，
$n = 0,1 \cdots ,N - 1$ ，$N$ 为样点数；$L$ 为分量数；${A_l}$ 、${f_l}$ 和${\varphi _l}$ 为第$l$ 个分量${s_l}\left[ n \right]$ 的复幅度、频率和初相，${f_s}$ 为采样率；$w\left[ n \right]$ 为加性高斯白噪声(additive white gaussian noise, AWGN)，$w\left[ n \right] \sim {\rm{CN}}( {0,2{\sigma ^2}} )$ 。为了获得接收信号频谱的更多细节，可以利用带补零的DFT变换，即：
$$R\left[ k \right] = \sum\limits_{n = 0}^{N - 1} {r\left[ n \right]{{\rm{e}}^{ - {\rm{j}}2{\text π} \frac{{kn}}{{N'}}}}} $$ (2) 式中，
$k = 0,1 \cdots ,N' - 1$ ；$N' = N + {P_0}$ ，${P_0}$ 为补零长度。将式(1)代入式(2)中，可得：
$$R\left[ k \right] = \sum\limits_{l = 1}^L {{S_l}\left[ k \right]} + \sum\limits_{n = 0}^{N - 1} {w\left[ n \right]{{\rm{e}}^{ - {\rm{j}}2{\text π} \frac{{kn}}{{N'}}}}} $$ (3) 式中，
$${S_l}\left[ k \right] = {A_l}{{\rm{e}}^{{\rm{j}}{\text π} \frac{{N - 1}}{{N'}}\left( {{k_l} + {\delta _l} - k} \right)}}\frac{{\sin \dfrac{{{\text π} N\left( {{k_l} + {\delta _l} - k} \right)}}{{N'}}}}{{\sin \dfrac{{{\text π} \left( {{k_l} + {\delta _l} - k} \right)}}{{N'}}}}$$ (4) $${k_l} = \left\lfloor {\left( {{{N'{f_l}} / {{f_s}}}} \right) + 0.5} \right\rfloor $$ (5) $${\delta_l} = \left( {{{N'{f_l}} / {{f_s}}}} \right) - {k_l}$$ (6) 根据式(4)可得：
$$S_l^{\left( 1 \right)}\sin \frac{{{\text π} {\delta _l} - {\text π} }}{{N'}} + S_l^{\left( { - 1} \right)}\sin \frac{{{\text π} {\delta _l} + {\text π} }}{{N'}} = \frac{{2S_l^{\left( 0 \right)}\sin \dfrac{{{\text π} {\delta _l}}}{{N'}}}}{{\sec \dfrac{{{\text π} N}}{{N'}}}}$$ (7) 式中，
$$S_l^{\left( x \right)} = {S_l}\left[ {{k_l} + x} \right]{{\rm{e}}^{{\rm{j}}{\text π} \frac{{N - 1}}{{N'}}x}}\text{,}\;\;x = - 1,0,1$$ (8) 通过分析式(7)，可得：
$${ \delta _l} = \frac{{N'}}{{\text π} }\arctan \frac{{\left( {S_l^{\left( { - 1} \right)} - S_l^{\left( 1 \right)}} \right)\tan \dfrac{{\text π} }{{N'}}}}{{2S_l^{\left( 0 \right)}\cos \dfrac{{{\text π} N}}{{N'}}\sec \dfrac{{\text π} }{{N'}} - S_l^{\left( { - 1} \right)} - S_l^{\left( 1 \right)}}}$$ (9) 当
$N' \gg \pi $ 时，式(9)可简化为：$${ \delta _l} \approx \Re \left\{ {\frac{{S_l^{\left( { - 1} \right)} - S_l^{\left( 1 \right)}}}{{2S_l^{\left( 0 \right)}\cos \dfrac{{{\text π} N}}{{N'}} - S_l^{\left( { - 1} \right)} - S_l^{\left( 1 \right)}}}} \right\}$$ (10) -
基于式(10)，本文提出一种加窗插值器，其推导过程如下。
为了简化分析，假设
$L = 1$ 。加窗的接收信号可以表示为：$$\tilde r\left[ n \right] = h\left[ n \right]r\left[ n \right]$$ (11) 式中，
$h\left[ n \right]$ 为窗函数。类似地，$\tilde r\left[ n \right]$ 的频域信号为：$$\tilde R\left[ k \right] = {\tilde S_l}\left[ k \right] + \sum\limits_{n = 0}^{N - 1} {h\left[ n \right]w\left[ n \right]{{\rm{e}}^{ - {\rm{j}}2{\text π} \frac{{kn}}{{N'}}}}} $$ (12) 式中，
$${\tilde S_l}\left[ k \right] = {A_l}\sum\limits_{n = 0}^{N - 1} {h\left[ n \right]{{\rm{e}}^{{\rm{j}}2{\text π} \frac{{{k_l} + {\delta _l} - k}}{{N'}}n}}} $$ (13) 根据式(10)，本文提出一种假设，即存在一个系数
${C_{N'}}$ ，使得式(14)为真。$${\delta _l} \approx {C_{N'}}\Re \left\{ {\frac{{\tilde S_l^{\left( { - 1} \right)} - \tilde S_l^{\left( 1 \right)}}}{{2\tilde S_l^{\left( 0 \right)}\cos \dfrac{{{\text π} N}}{{N'}} - \tilde S_l^{\left( { - 1} \right)} - \tilde S_l^{\left( 1 \right)}}}} \right\}$$ (14) 式中，
$$\tilde S_l^{\left( x \right)} = {\tilde S_l}\left[ {{k_l} + x} \right]{{\rm{e}}^{{\rm{j}}{\text π} \frac{{N - 1}}{{N'}}x}},x = - 1,0,1$$ (15) 根据式(13)可知：
$${\tilde S_l}\left[ {{k_l} - x} \right] = {A_l}f\left( {x + { \delta_l}} \right)$$ (16) 式中，
$$f\left( x \right) = \sum\limits_{n = 0}^{N - 1} {h\left[ n \right]{{\rm{e}}^{{\rm{j}}2{\text π} \frac{{xn}}{{N'}}}}} $$ (17) 根据泰勒展开可知，
$f\left( {x + {\delta _l}} \right)$ 可展开为：$$f\left( {x + {\delta_l}} \right) = f\left( x \right) + f'\left( x \right){ \delta _l} + O\left( {\delta _l^2} \right)$$ (18) 式中，
$$f'\left( x \right) = \frac{{{\rm{j}}2{\text π} }}{{N'}}\sum\limits_{n = 0}^{N - 1} {nh\left[ n \right]{{\rm{e}}^{{\rm{j}}2{\text π} \frac{{xn}}{{N'}}}}} $$ (19) 根据式(15)、式(16)和式(18)可知，式(14)中的分子可以展开为：
$$\tilde S_l^{\left( { - 1} \right)} - \tilde S_l^{\left( 1 \right)} = {A_l}\left[ {{\rm{j}}{a_0} + {a_1}{\delta _l} + O\left( {\delta _l^2} \right)} \right]$$ (20) 式中，
$${a_0} = \Im \left\{ {f\left( 1 \right){{\rm{e}}^{ - {\rm{j}}{\text π} \frac{{N - 1}}{{N'}}}} - f\left( { - 1} \right){{\rm{e}}^{{\rm{j}}{\text π} \frac{{N - 1}}{{N'}}}}} \right\}$$ (21) $${a_1} = \Re \left\{ {f'\left( 1 \right){{\rm{e}}^{ - {\rm{j}}{\text π} \frac{{N - 1}}{{N'}}}} - f'\left( { - 1} \right){{\rm{e}}^{{\rm{j}}{\text π} \frac{{N - 1}}{{N'}}}}} \right\}$$ (22) 类似地，式(14)中的分母可以展开为：
$$\frac{{2\tilde S_l^{\left( 0 \right)}}}{{\sec \dfrac{{{\text π} N}}{{N'}}}} - \tilde S_l^{\left( { - 1} \right)} - \tilde S_l^{\left( 1 \right)} = {A_l}\left[ {{b_0} + {\rm{j}}{b_1}{ \delta _l} + O\left( { \delta _l^2} \right)} \right]$$ (23) 式中，
$${b_0} = \Re \left\{ {2f\left( 0 \right)\cos \frac{{{\text π} N}}{{N'}} - \frac{{f\left( 1 \right)}}{{{{\rm{e}}^{{\rm{j}}{\text π} \frac{{N - 1}}{{N'}}}}}} - \frac{{f\left( { - 1} \right)}}{{{{\rm{e}}^{ - {\rm{j}}{\text π} \frac{{N - 1}}{{N'}}}}}}} \right\}$$ (24) $${b_1} = \Im \left\{ {2f'\left( 0 \right)\cos \frac{{{\text π} N}}{{N'}} - \frac{{f'\left( 1 \right)}}{{{{\rm{e}}^{{\rm{j}}{\text π} \frac{{N - 1}}{{N'}}}}}} - \frac{{f'\left( { - 1} \right)}}{{{{\rm{e}}^{ - {\rm{j}}{\text π} \frac{{N - 1}}{{N'}}}}}}} \right\}$$ (25) 因为
${a_0}$ 、${a_1}$ 、${b_0}$ 和${b_1}$ 都是实数，所以有：$$\Re \left\{ {\frac{{{\rm{j}}{a_0} + {a_1}{ \delta_l}}}{{{b_0} + {\rm{j}}{b_1}{ \delta_l}}}} \right\} = \frac{{{a_1}{b_0} + {a_0}{b_1}}}{{b_0^2}}{\delta_l} + O\left( {\delta_l^2} \right)$$ (26) 将式(20)和式(23)代入式(26)可得：
$${ \delta _l} = \Re \left\{ {\frac{{\dfrac{{b_0^2}}{{{a_1}{b_0} + {a_0}{b_1}}}\left( {\tilde S_l^{\left( { - 1} \right)} - \tilde S_l^{\left( 1 \right)}} \right)}}{{2\tilde S_l^{\left( 0 \right)}\cos \dfrac{{ {\text π} N}}{{N'}} - \tilde S_l^{\left( { - 1} \right)} - \tilde S_l^{\left( 1 \right)}}}} \right\} + O\left( {\delta_l^2} \right)$$ (27) 根据式(27)可知，本小节所提出的假设是真，且系数
${C_{N'}}$ 可以表示为：$${C_{N'}} = \frac{{b_0^2}}{{{a_1}{b_0} + {a_0}{b_1}}}$$ (28) 特别地，当窗函数为矩形窗时，有
${C_{N'}} = 1$ 。基于式(27)，本文提出了一种加窗插值器，即：
$${\hat \delta _l} = {C_{N'}}\Re \left\{ {\frac{{\tilde R_l^{\left( { - 1} \right)} - \tilde R_l^{\left( 1 \right)}}}{{2\tilde R_l^{\left( 0 \right)}\cos \dfrac{{ {\text π} N}}{{N'}} - \tilde R_l^{\left( { - 1} \right)} - \tilde R_l^{\left( 1 \right)}}}} \right\}$$ (29) 式中，
${\hat \delta _l}$ 为${\delta _l}$ 的估计值，且：$$\tilde R_l^{\left( x \right)} = \tilde R\left[ {{{\hat k}_l} + x} \right]{{\rm{e}}^{{\rm{j}} {\text π} \frac{{N - 1}}{{N'}}x}},x = - 1,0,1$$ (30) 式中，
${\hat k_l}$ 为${k_l}$ 的估计值，具体运算在下一节给出。 -
不失一般性地，假设
$\left| {{A_1}} \right| \geqslant \left| {{A_2}} \right| \geqslant \cdots \geqslant \left| {{A_L}} \right|$ 。在这个阶段，$l$ 被依次设置为$1,2, \cdots ,L$ ，当前$l - 1$ 个分量被估计后，可以根据式(31)将这些分量从接收信号中消除以降低分量间互干扰：$${r_l}\left[ n \right] = r\left[ n \right] - \sum\limits_{l' = 1}^{l - 1} {{{\hat A}_{l'}}{{\rm{e}}^{{\rm{j}}2 {\text π} \frac{{{{\hat f}_{l'}}}}{{{f_s}}}n}}} $$ (31) 根据文献[16]可知，
${k_l}$ 可估计为：$${\hat k_l} = \arg \mathop {\max }\limits_{0 \leqslant k < N'} \left\{ {{{\left| {{R_l}\left[ k \right]} \right|}^2}} \right\}$$ (32) 式中，
$${R_l}\left[ k \right] = \sum\limits_{n = 0}^{N - 1} {{r_l}\left[ n \right]{{\rm{e}}^{ - {\rm{j}}2 {\text π} \frac{{kn}}{{N'}}}}} $$ (33) 因为
${r_l}\left[ n \right]$ 等同于被加了一个矩形窗，所以可将$\tilde R\left[ k \right] = {R_l}\left[ k \right]$ 代入式(30)得到$\tilde R_l^{\left( x \right)}$ ，然后将$\tilde R_l^{\left( x \right)}$ 和${C_{N'}} = $ $ 1$ 代入式(29)得到${\hat \delta _l}$ ，这样，${f_l}$ 可估计为：$${\hat f_{l,1}} = \frac{{{{\hat k}_l} + {{\hat \delta }_l}}}{{N'}}{f_s}$$ (34) 为了降低分量间互干扰和高阶项
$O( {\delta _l^2} )$ 对频率估计精度的影响，可以通过${\hat f_{l,1}}$ 将第$l$ 个分量的频率移至零附近以减小高阶项，并对接收信号采用能降低频谱泄露的非矩形窗，即：$${{\overset{\smile} r} _l}\left[ n \right] = h\left[ n \right]{r_l}\left[ n \right]{{\rm{e}}^{ - {\rm{j}}2 {\text π} \frac{{{{\hat f}_{l,1}}}}{{{f_s}}}n}}$$ (35) 根据上一节的分析可知，
${f_l}$ 可被重新估计为：$${\hat f_l} = {\hat f_{l,1}} + \Delta {\hat f_l}$$ (36) 式中，
$\Delta {\hat f_l}$ 为$\Delta {f_l} = {f_l} - {\hat f_{l,1}}$ 的估计值，且：$$\Delta {\hat f_l} = \frac{{{f_s}}}{{N'}}{C_{N'}}\Re \left\{ {\frac{{{\overset{\smile} R} _l^{\left( { - 1} \right)} - {\overset{\smile} R} _l^{\left( 1 \right)}}}{{2{\overset{\smile} R} _l^{\left( 0 \right)}\cos \dfrac{{{\text π} N}}{{N'}} - {\overset{\smile} R} _l^{\left( { - 1} \right)} - {\overset{\smile} R} _l^{\left( 1 \right)}}}} \right\}$$ (37) $${\overset{\smile} R} _l^{\left( x \right)} = {{\rm{e}}^{{\rm{j}}{\text π} \frac{{N - 1}}{{N'}}x}}\sum\limits_{n = 0}^{N - 1} {{\overset{\smile} r} \left[ n \right]{{\rm{e}}^{ - {\rm{j}}2{\text π} \frac{{xn}}{{N'}}}}}\text{,}\;\;x = - 1,0,1$$ (38) 此外，
${A_l}$ 可估计为：$${\hat A_l} = \frac{1}{N}\sum\limits_{n = 0}^{N - 1} {{r_l}\left[ n \right]{{\rm{e}}^{ - {\rm{j}}2{\text π} \frac{{{{\hat f}_l}}}{{{f_s}}}n}}} $$ (39) 在这一阶段的最后，可得到所有参数的估计值。
-
虽然第一阶段的非矩形窗可以降低频谱泄露，但是也会造成一定的SNR损失，因此第一阶段的频率估计值的精度可以被进一步提升。
在这一阶段中，
$l$ 仍然被依次设置为$1,2, \cdots ,L$ 。因为在第一阶段所有分量的参数估计值都已得到，所以分量间的互干扰可以被消除为：$${\dot r_l}\left[ n \right] = r\left[ n \right] - \sum\limits_{l' = 1,l' \ne l}^L {{{\hat A}_{l'}}{{\rm{e}}^{{\rm{j}}2{\text π} \frac{{{{\hat f}_{l'}}}}{{{f_s}}}n}}} $$ (40) 第二阶段采用矩形窗，并可将
${C_{N'}} = 1$ 、$h\left[ n \right] = $ $ 1$ 、${\hat f_{l,1}} = {\hat f_l}$ 和${r_l}\left[ n \right] = {\dot r_l}\left[ n \right]$ 代入式(35)～式(38)更新${\hat f_l}$ ，然后将${r_l}\left[ n \right] = {\dot r_l}\left[ n \right]$ 和${\hat f_l}$ 代入式(39)更新${\hat A_l}$ 。在第二阶段的最后，所有参数估计值都可更新。基于TSWI的多音信号频率估计算法的详细步骤如下一节所示。注意，为了降低计算复杂度，补零长度
${P_0}$ 首先设置为${2^{\left\lceil {{{\log }_2}N} \right\rceil }} - N$ ，此时DFT可以通过快速傅里叶变换(fast Fourier transform, FFT)实现，然后为了提高估计精度，${P_0}$ 设置为$M$ ，$M$ 为一个远大于$N$ 的整数。 -
根据基于TSWI的多音信号频率估计算法可知，频率估计性能主要由
$\Delta {\hat f_l}$ 决定。输入：
$r\left[ n \right]$ 、非矩形窗$h\left[ n \right]$ 及其对应的系数${C_{N + M}}$ 输出：
${\hat f_l}$ ，其中$l = 1,2, \cdots ,L$ 第一阶段：
for
$l = 1,2, \cdots ,L$ 根据式(31)消除
$r\left[ n \right]$ 中前$l - 1$ 个分量，得到${r_l}\left[ n \right]$ ；在
${r_l}\left[ n \right]$ 后补${2^{\left\lceil {{{\log }_2}N} \right\rceil }} - N$ 个零，计算$N' = {2^{\left\lceil {{{\log }_2}N} \right\rceil }}$ ；对补零后的
${r_l}\left[ n \right]$ 进行FFT得到${R_l}\left[ k \right]$ ；将
${R_l}\left[ k \right]$ 代入式(32)得到${\hat k_l}$ ；将
$\tilde R\left[ k \right] = {R_l}\left[ k \right]$ 代入式(30)得到$\tilde R_l^{\left( x \right)}$ ，其中x =$ \pm 1,0$ ；将
$\tilde R_l^{\left( x \right)}$ 和${C_{N'}} = 1$ 代入式(29)得到${\hat \delta _l}$ ；将
${\hat k_l}$ 和${\hat \delta _l}$ 代入式(34)得到${\hat f_{l,1}}$ ；更新
$N' = N + M$ ，${C_{N'}} = {C_{N + M}}$ ；将
${\hat f_{l,1}}$ 、${r_l}\left[ n \right]$ 和$h\left[ n \right]$ 代入式(35)得到${{\overset{\smile} r} _l}\left[ n \right]$ ；将
${{\overset{\smile} r} _l}\left[ n \right]$ 代入式(38)得到${\overset{\smile} R} _l^{\left( x \right)}$ ，其中$ x= \pm 1,0$ ；将
${\overset{\smile} R} _l^{\left( x \right)}$ 、${C_{N'}}$ 和${\hat f_{l,1}}$ 代入式(37)得到$\Delta {\hat f_l}$ ；将
${\hat f_{l,1}}$ 和$\Delta {\hat f_l}$ 代入式(36)得到${\hat f_l}$ ；将
${r_l}\left[ n \right]$ 和${\hat f_l}$ 代入式(39)得到${\hat A_l}$ ；end
第二阶段：
更新
$h\left[ n \right] = 1$ ，${C_{N'}} = 1$ ；for
$l = 1,2, \cdots ,L$ 根据式(40)消除其他分量，得到
${\dot r_l}\left[ n \right]$ ；将
${r_l}\left[ n \right] = {\dot r_l}\left[ n \right]$ 和${\hat f_{l,1}} = {f_l}$ 代入式(35)得到${{\overset{\smile} r} _l}\left[ n \right]$ ；将
${{\overset{\smile} r} _l}\left[ n \right]$ 代入式(38)得到${\overset{\smile} R} _l^{\left( x \right)}$ ，其中$x = \pm 1,0$ ；将
${\overset{\smile} R} _l^{\left( x \right)}$ 和${\hat f_{l,1}} = {f_l}$ 代入式(37)得到$\Delta {\hat f_l}$ ；将
${\hat f_{l,1}} = {f_l}$ 和$\Delta {\hat f_l}$ 代入式(36)更新${\hat f_l}$ ；将
${r_l}\left[ n \right]$ 和${\hat f_l}$ 代入式(39)更新${\hat A_l}$ ；end
根据式(37)可知：
$$\Delta {\hat f_l} \approx \frac{{{f_s}}}{{N'}}{C_{N'}}\Re \left\{ {\frac{{{\overset{\smile} S} _l^{\left( { - 1} \right)} - {\overset{\smile} S} _l^{\left( { - 1} \right)} + {W_1}}}{{S + {W_2}}}} \right\}$$ (41) 式中，
$$S = 2{\overset{\smile} S} _l^{\left( 0 \right)}\cos \frac{{{\text π} N}}{{N'}} - {\overset{\smile} S} _l^{\left( { - 1} \right)} - {\overset{\smile} S} _l^{\left( 1 \right)}$$ (42) $${\overset{\smile} S} _l^{\left( x \right)} = {A_l}{{\rm{e}}^{{\rm{j}}{\text π} \frac{{N - 1}}{{N'}}x}}\sum\limits_{n = 0}^{N - 1} {h\left[ n \right]{{\rm{e}}^{{\rm{j}}2{\text π} \left( {\frac{{\Delta {f_l}}}{{{f_s}}} - \frac{x}{{N'}}} \right)n}}} $$ (43) $${W_1} = 2j\sum\limits_{n = 0}^{N - 1} {{\overset{\smile} w} \left[ n \right]\sin \frac{{2{\text π} ( {n - \tilde N} )}}{{N'}}} $$ (44) $${W_2} = 2\sum\limits_{n = 0}^{N - 1} {{\overset{\smile} w} \left[ n \right]\left[ {\cos \frac{{{\text π} N}}{{N'}} - \cos \frac{{2{\text π} ( {n - \tilde N} )}}{{N'}}} \right]} $$ (45) $${\overset{\smile} w} \left[ n \right] = h\left[ n \right]w\left[ n \right]{{\rm{e}}^{ - {\rm{j}}2{\text π} \frac{{{{\hat f}_{l,1}}}}{{{f_s}}}n}}$$ (46) $$\tilde N = {{\left( {N - 1} \right)} / 2}$$ (47) 当SNR较高时，认为
${{\Delta {f_l}} / {{f_s}}} \approx 0$ ，$S$ 简化为：$$S = 2{A_l}\sum\limits_{n = 0}^{N - 1} {h\left[ n \right]\left[ {\cos \frac{{{\text π} N}}{{N'}} - \cos \frac{{2{\text π} ( {n - \tilde N} )}}{{N'}}} \right]} $$ (48) 同理，式(41)可以简化为：
$$\Delta {\hat f_l} = \Delta {f_l} + \frac{{\Re \left\{ {\dfrac{{j{f_s}{C_{N'}}}}{{{A_l}N'}}\displaystyle\sum\limits_{n = 0}^{N - 1} {{\overset{\smile} w} \left[ n \right]\sin \dfrac{{2{\text π} ( {n - \tilde N} )}}{{N'}}} } \right\}}}{{\displaystyle\sum\limits_{n = 0}^{N - 1} {h\left[ n \right]\left[ {\cos \dfrac{{{\text π} N}}{{N'}} - \cos \dfrac{{2{\text π} ( {n - \tilde N} )}}{{N'}}} \right]} }}$$ (49) 因为有
$\Delta {f_l} = {f_l} - {\hat f_{l,1}}$ 且${\hat f_l} = {\hat f_{l,1}} + \Delta {\hat f_l}$ ，所以有：$${\hat f_l} - {f_l} = \Delta {\hat f_l} - \Delta {f_l}$$ (50) 将式(49)代入式(50)可得：
$${\hat f_l} = {f_l} + \frac{{\Re \left\{ {\dfrac{{{\rm{j}}{f_s}{C_{N'}}}}{{{A_l}N'}}\displaystyle\sum\limits_{n = 0}^{N - 1} {{\overset{\smile} w} \left[ n \right]\sin \dfrac{{2{\text π} ( {n - \tilde N} )}}{{N'}}} } \right\}}}{{\displaystyle\sum\limits_{n = 0}^{N - 1} {h\left[ n \right]\left[ {\cos \dfrac{{{\text π} N}}{{N'}} - \cos \frac{{2{\text π} ( {n - \tilde N} )}}{{N'}}} \right]} }}$$ (51) 根据式(46)可知：
$${\overset{\smile} w} \left[ n \right] \sim {\rm{CN}}( {0,2{h^2}\left[ n \right]{\sigma ^2}} )$$ (52) 根据式(51)和式(52)可知，
${\hat f_l}$ 的均方误差(mean-square error, MSE)可表示为:$$\sigma _l^2 \approx \frac{{\dfrac{{f_s^2C_{N'}^2}}{{2{\rm{SN}}{{\rm{R}}_l}{{N'}^2}}}\displaystyle\sum\limits_{n = 0}^{N - 1} {{h^2}\left[ n \right]{{\sin }^2}\dfrac{{2{\text π} ( {n - \tilde N} )}}{{N'}}} }}{{{{\left\{ {\displaystyle\sum\limits_{n = 0}^{N - 1} {h\left[ n \right]\left[ {\cos \dfrac{{{\text π} N}}{{N'}} - \cos \dfrac{{2{\text π} ( {n - \tilde N} )}}{{N'}}} \right]} } \right\}}^2}}}$$ (53) 式中，
${\rm{SN}}{{\rm{R}}_l} = {{{{\left| {{A_l}} \right|}^2}} / {\left( {2{\sigma ^2}} \right)}}$ 。由于第二阶段采用矩形窗，式(53)可以简化为：
$$\sigma _l^2 \approx \frac{{\frac{{f_s^2C_{N'}^2}}{{2{{N'}^2}}}\left( {\dfrac{N}{2} - \dfrac{{N'}}{{4{\text π} }}\sin \dfrac{{2{\text π} N}}{{N'}}} \right)}}{{{\rm{SN}}{{\rm{R}}_l}{{\left( {N\cos \dfrac{{{\text π} N}}{{N'}} - \dfrac{{N'}}{{\text π} }\sin \dfrac{{{\text π} N}}{{N'}}} \right)}^2}}}$$ (54) -
本节仿真了不同多音信号频率估计算法，并对其性能进行了对比分析，仿真参数如下：
$L = 2$ ，$\left| {{A_1}} \right| = \left| {{A_2}} \right| = 1$ ，${f_1} \sim {\rm{U}}\left[ {0,{f_s}} \right)$ ，${f_2} = {f_1} + \alpha $ ，$\alpha $ 为频率间隔，$\angle {A_1} \sim {\rm{U}}$ $ \left[ { -{\text {π}} ,{\text π} } \right) $ ，$\angle {A_2} \sim {\rm{U}}$ $ \left[ { -{\text {π}} ,{\text π} } \right) $ ，$M = 3N$ ，$N = $ $ 256$ ，沿用文献[13]算法的迭代次数为2。图1给出了TSWI第一阶段采用不同窗函数时的两个分量的平均均方MSE (root MSE，RMSE)性能，其中
$\alpha = {{2.4{f_s}} / N}$ ，一个DFT栅格为${{{f_s}} / N}$ ，${\rm{SNR}}$ 为所有分量的信噪比，其与${\rm{SN}}{{\rm{R}}_l}$ 的关系可以表示为：$${\rm{SN}}{{\rm{R}}_l} = {{{{\left| {{A_l}} \right|}^2}{\rm{SNR}}} / {\sum\nolimits_{l'}^L {{{\left| {{A_{l'}}} \right|}^2}} }}$$ (55) 从图1中可以看出，采用非矩形窗时的频率估计性能优于采用矩形窗时的频率估计性能，这是因为这些非矩形窗可以降低频率泄露，减少分量间的互干扰。此外，采用汉明窗时TSWI具有最优的频率估计性能，当
${\rm{SNR}} \geqslant - 5\;{\rm{dB}}$ 时，其RMSE仿真曲线与理论RMSE曲线基本重合，且与CRLB也基本重合，在之后的仿真中，TSWI第一阶段都采用汉明窗。
图 1 TSWI第一阶段采用不同窗函数时的两分量平均RMSE
图2展示了
$\alpha = {{2.4{f_s}} / N}$ 时不同多音信号频率估计算法的平均RMSE，由于频率间隔较小，分量间互干扰较大，文献[11]算法、文献[13]算法和抛物线插值算法都出现估计误差平层。与这些算法不同，TSWI未出现估计误差平层，当${\rm{SNR}} \geqslant - 5\;{\rm{dB}}$ 时，其RMSE基本达到了CRLB，这是因为TSWI在第一阶段通过汉明窗抑制了频谱泄露，并在第二阶段通过矩形窗弥补了第一阶段汉明窗造成的SNR损失。
图 2
$\alpha = {{2.4{f_s}} / N}$ 时不同多音频率估计算法的平均RMSE在图3中，
$\alpha $ 被设置为${{31.9{f_s}} / N}$ ，由于频率间隔较大，所有算法都未出现频率估计误差平层，然而，TSWI的估计性能仍然优于其他算法。
图 3
$\alpha = {{31.9{f_s}} / N}$ 时不同多音频率估计算法的平均RMSE图4对比了不同频率间隔下各多音信号频率估计算法的平均RMSE，其中SNR=30 dB。可以发现，TSWI具有最优的频率估计性能，当
$\alpha > {{2{f_s}} / N}$ 时，TSWI的RMSE基本达到了CRLB，其他算法都需要更大的频率间隔才能不再受分量间互干扰的影响。
图 4 SNR=30 dB时不同多音频率估计算法的平均RMSE
图5将
$L$ 增大到4，各分量的频率分别为${{5.1{f_s}} / N}$ 、${{7.2{f_s}} / N}$ 、${{16.3{f_s}} / N}$ 和${{35.8{f_s}} / N}$ ，各分量的幅度都为1，且初相都服从区间为$ -{\text π} $ ～${\text π} $ 的均匀分布。由于$L$ 增大到4，所有算法都出现了频率估计平层，然而TSWI在SNR>30 dB时才出现误差平层，其他算法在SNR>20 dB时就出现了误差平层，此外，TSWI的误差平层远小于其他算法。
图 5 L=4时不同多音频率估计算法的平均RMSE
-
本文提出了一种基于TSWI的多音信号频率估计算法，本算法包含一个新的插值器，其插值器支持任意窗函数，通过在第一阶段选择可以抑制频谱泄露的非矩形窗并在第二阶段选择矩形窗，在降低分量间互干扰的同时不损失SNR。仿真结果表明，本文算法具有比现有算法更优的频率估计性能。
Frequency Estimation of Multi-Tone by Two-Stage Windowed Interpolation
-
摘要: 频域插值是一种广泛应用于多音信号频率估计的方法。为了提高相邻单音分量频率间隔较小时的频率估计性能，该文提出了一种基于两阶段加窗插值的频率估计算法。该算法采用一种新的支持任意窗函数的插值器来估计频率，通过在不同的阶段选择不同的窗函数，可以在不损失信噪比的前提下减少多个单音分量之间的相互干扰。数值结果表明，该算法具有比现有算法更好的估计性能，特别是在相邻单音分量频率间隔较小的情况下。Abstract: Frequency-domain interpolation is a widely used method for frequency estimation of multiple sinusoids. To improve the frequency estimation performance when the interval between adjacent frequencies is small, a new estimation algorithm based on two-stage windowed interpolation is proposed in this paper. In this algorithm, frequencies are estimated by a new interpolator which supports arbitrary windows. By choosing different windows in different stages, the proposed algorithm can reduce the mutual interference between multiple sinusoids without loss of signal-to-noise ratio. Numerical results show that the proposed frequency estimation algorithm has better estimation performance than the state-of-the-art algorithms, especially when the interval between adjacent frequencies is small.
-
Key words:
- discrete Fourier transforms /
- frequency estimation /
- interpolation /
- windows
-
[1] DING Y P, SUN Y H, HUANG G W, et al. Human target localization using Doppler through-wall radar based on micro-Doppler frequency estimation[J]. IEEE Sensors Journal, 2020, 20(15): 8778-8788. doi: 10.1109/JSEN.2020.2983104 [2] 高尚蕾, 张治中, 段浴, 等. 5G系统基于PSS和SSS联合频偏估计方法[J]. 光通信研究, 2020(6): 65-69. GAO S L, ZHANG Z Z, DUAN Y, et al. A joint frequency offset estimation method for 5G system based on PSS and SSS[J]. Study on Optical Communications, 2020(6): 65-69. [3] 何金埔. NC-OFDM系统干扰对消技术研究与实现[D]. 成都: 电子科技大学, 2020. HE J P. Research and implementation of jamming cancellation technology in NC-OFDM system[D]. Chengdu: University of Electronic Science and Technology of China, 2020. [4] IVANOV D, ZHDANOV A. Weighted total least squares for frequency estimation of real sinusoids based on augmented system[C]//2020 IEEE East-West Design & Test Symposium. Varna, Bulgaria: IEEE, 2020: 1-5. [5] SERBES A. Fast and efficient sinusoidal frequency estimation by using the DFT coefficients[J]. IEEE Trans on Communications, 2019, 67(3): 2333-2342. doi: 10.1109/TCOMM.2018.2886355 [6] BELEGA D, PETRI D, DALLET D. Frequency estimation of a sinusoidal signal via a three-point interpolated DFT method with high image component interference rejection capability[J]. Digital Signal Processing, 2014, 24: 162-169. doi: 10.1016/j.dsp.2013.09.014 [7] DIAO R P, MENG Q F. An interpolation algorithm for discrete Fourier transforms of weighted damped sinusoidal signals[J]. IEEE Trans on Instrumentation and Measurement, 2014, 63(6): 1505-1513. doi: 10.1109/TIM.2013.2289585 [8] BELEGA D, PETRI D. Frequency estimation by two- or three-point interpolated Fourier algorithms based on cosine windows[J]. Signal Processing, 2015, 117: 115-125. doi: 10.1016/j.sigpro.2015.05.005 [9] LUO J F, XIE Z J, XIE M. Interpolated DFT algorithms with zero padding for classic windows[J]. Mechanical Systems and Signal Processing, 2016, 70: 1011-1025. [10] DUDA K. DFT interpolation algorithm for Kaiser-Bessel and Dolph-Chebyshev windows[J]. IEEE Trans on Instrumentation and Measurement, 2011, 60(3): 784-790. doi: 10.1109/TIM.2010.2046594 [11] CANDAN C. Fine resolution frequency estimation from three DFT samples: Case of windowed data[J]. Signal Processing, 2015, 114: 245-250. doi: 10.1016/j.sigpro.2015.03.009 [12] MAMANDIPOOR B, RAMASAMY D, MADHOW U. Newtonized orthogonal matching pursuit: Frequency estimation over the continuum[J]. IEEE Trans on Signal Processing, 2016, 64(19): 5066-5081. doi: 10.1109/TSP.2016.2580523 [13] YE S L, ABOUTANIOS E. Rapid accurate frequency estimation of multiple resolved exponentials in noise[J]. Signal Processing, 2017, 132: 29-39. doi: 10.1016/j.sigpro.2016.09.010 [14] ABOUTANIOS E, MULGREW B. Iterative frequency estimation by interpolation on Fourier coefficients[J]. IEEE Trans on Signal Processing, 2005, 53(4): 1237-1242. doi: 10.1109/TSP.2005.843719 [15] DJUKANOVIC S, POPOVIC V. Efficient and accurate detection and frequency estimation of multiple sinusoids[J]. IEEE Access, 2019, 7: 1118-1125. doi: 10.1109/ACCESS.2018.2886397 [16] DJUKANOVIC S, POPOVIC T, MITROVIC A. Precise sinusoid frequency estimation based on parabolic interpolation[C]//Telecommunications Forum. Belgrade, Serbia: IEEE, 2016: 1-4. -
点击查看大图
图(5)
计量
- 文章访问数: 4635
- HTML全文浏览量: 1527
- PDF下载量: 41
- 被引次数: 0
