CN106533394A - 一种基于自适应滤波器幅频响应的高精度频率估计方法 - Google Patents

Abstract

一种基于自适应滤波器幅频响应的高精度频率估计方法，包括计算信号协方差矩阵、确定第一个正弦信号的频率区间、计算自适应滤波器的权系数向量、计算自适应滤波器的幅频响应、计算信号伪功率谱、计算第一个正弦信号的频率估计值。本发明在同样的信噪比及信号采样点数等条件下，本发明方法的频率估计精度优于或接近已有频率估计方法的频率估计精度，特别是在信噪比很低、信号采样点数很少等情况下，相比于已有的频率估计方法，本发明方法能够在正弦信号频率附近得到更加明显的谱峰，从而得到更加精确的频率估计值。

Description

一种基于自适应滤波器幅频响应的高精度频率估计方法

技术领域

本发明涉及信号谱估计领域，更具体地说，涉及一种基于自适应滤波器幅频响应的高精度频率估计方法。

背景技术

信号谱估计是指对一个有限长的平稳序列，估计其在整个频率域内的功率分布。信号谱估计在通信、雷达、声呐、地震学等诸多领域中均有广泛应用。以雷达和声呐系统为例，人们可以根据接收信号的频谱对观测范围内的目标进行定位。根据功率谱的分布情况，信号谱估计可以分为两类即：连续谱估计和线谱估计，其中线谱估计的主要任务就是对包含正弦信号和噪声的观测信号进行处理，以获得正弦信号的频率估计值。

已有的频率估计方法包括：快速傅里叶变换(Fast Fourier Transformation,FFT)方法、多重信号分类(Multiple Signal Classification,MUSIC)方法、求根多重信号分类(Root-MUSIC)方法、信号参数估计旋转不变技术(Estimating Signal Parametersvia Rotational Invariance Techniques,ESPRIT)等。其中，FFT方法计算量小且易于实现，因此在实际工程中得到了广泛的应用，但是FFT方法的频率分辨率不够高；MUSIC方法、Root-MUSIC方法以及ESPRIT方法属于子空间处理方法，它们在高信噪比情况下可获得更高的频率分辨率。但是，当信噪比很低且信号采样点数很少时，不管是常规的FFT方法，还是MUSIC、Root-MUSIC及ESPRIT等子空间方法，其估计性能均显著恶化，所得到的频率估计精度远远不能达到实际应用的要求。因此，研究一种适用于信噪比很低、信号采样点数很少等情况下的高精度频率估计方法具有重要的实用价值。

发明内容

为了有效解决已有频率估计方法在信噪比很低且信号采样点数很少等情况下频率估计性能显著恶化的问题，本发明提供了一种基于自适应滤波器幅频响应的高精度频率估计方法。在信噪比很低和信号采样点数很少等情况下，本发明方法的频率估计精度明显优于已有频率估计方法。

本发明的基本思路是：首先，对观测的连续时间信号进行采样得到离散时间信号，并由离散时间信号计算得到信号协方差矩阵；其次，由先验信息确定第一个正弦信号的频率区间；再次，由信号协方差矩阵和频率区间计算自适应滤波器的权系数向量，并使自适应滤波器的幅频响应满足以下特点：当自适应滤波器的中心频率不等于第一个正弦信号的频率时，自适应滤波器的幅频响应在频率区间内存在明显的零陷，而当自适应滤波器的中心频率等于第一个正弦信号的频率时，自适应滤波器的幅频响应在频率区间内不存在明显的零陷；然后，由自适应滤波器的权系数向量计算自适应滤波器在频率区间内的幅频响应；接着，在频率区间内搜索自适应滤波器的幅频响应的最小值及最大值，并由上述最小值和最大值计算信号伪功率谱；最后，在频率区间内搜索信号伪功率谱的谱峰，并将信号伪功率谱的谱峰对应的频率作为第一个正弦信号的频率估计值。若在观测信号中正弦信号的个数大于1，则根据先验信息确定其它正弦信号的频率区间，并重复上述步骤即可得到其它正弦信号的频率估计值。

本发明的技术方案是：一种基于自适应滤波器幅频响应的高精度频率估计方法，具体包括下述步骤：

第一步：计算信号协方差矩阵

首先，对观测的连续时间信号进行采样得到离散时间信号{x(1),x(2),...,x(N)}，其中采样频率为f_s，单位为Hz，N表示离散时间信号的采样点数。然后，由离散时间信号{x(1),x(2),...,x(N)}构造信号向量x(k)，即：

x(k)＝[x(k-m+1),x(k-m+2),...,x(k)]^H

其中，m为信号向量x(k)的长度，m为正整数，且1＜m＜N，k也为正整数，且m≤k≤N，m的具体大小根据实际应用情况确定，符号[·]^H表示共轭转置操作。最后，由信号向量x(k)计算信号协方差矩阵R，即：

其中x^H(k)表示对信号向量x(k)的共轭转置。

第二步：确定第一个正弦信号的频率区间

根据先验信息确定第一个正弦信号的频率区间，并记为Φ₁＝[f₁ ^l,f₁ ^r]，其中f₁ ^l为频率区间的下边界，f₁ ^r为频率区间的上边界，单位均为Hz。第一个正弦信号的频率落在频率区间Φ₁之中。

第三步：计算自适应滤波器的权系数向量

根据第一步得到的信号协方差矩阵R以及第二步得到的第一个正弦信号的频率区间Φ₁，按照下式计算自适应滤波器的权系数向量h，即：

其中，R^-1表示信号协方差矩阵R的逆，向量f_c为自适应滤波器的中心频率，f_c∈Φ₁表示f_c的取值范围在频率区间Φ₁内。按照上述方式设计的自适应滤波器具有以下特性：当自适应滤波器的中心频率不等于第一个正弦信号的频率时，自适应滤波器的幅频响应在频率区间内存在明显的零陷，而当自适应滤波器的中心频率等于第一个正弦信号的频率时，自适应滤波器的幅频响应在频率区间内不存在明显的零陷。

第四步：计算自适应滤波器的幅频响应

根据第三步计算得到的自适应滤波器的权系数向量h，计算自适应滤波器在频率区间Φ₁内的幅频响应H(f)，即：

H(f)＝|a^H(f)h|,f∈Φ₁

其中，向量|·|表示求复数的幅度。自适应滤波器在频率区间Φ₁内的幅频响应H(f)也可以由信号协方差矩阵R和自适应滤波器的中心频率f_c直接计算得到，即：

第五步：计算信号伪功率谱

首先，根据第四步计算得到的自适应滤波器在频率区间Φ₁内的幅频响应H(f)，搜索其最大值和最小值，分别记为：

然后，由自适应滤波器幅频响应最大值H_max(f_c)和自适应滤波器幅频响应最小值H_min(f_c)计算信号伪功率谱P(f_c)，即：

信号伪功率谱P(f_c)也可以由信号协方差矩阵R和自适应滤波器的中心频率f_c直接计算得到，即：

第六步：计算第一个正弦信号的频率估计值

根据第五步计算得到的信号伪功率谱P(f_c)，搜索其在频率区间Φ₁内的谱峰，并将信号伪功率谱P(f_c)的谱峰对应的频率作为第一个正弦信号的频率估计值f₁，即：

若在观测信号中正弦信号的个数大于1，则首先根据先验信息确定其它正弦信号的频率区间，然后重复上述步骤中的第三步至第六步，即可得到其它正弦信号的频率估计值。

本发明的有益效果是：在同样的信噪比及信号采样点数等条件下，本发明方法的频率估计精度优于或接近已有频率估计方法的频率估计精度，特别是在信噪比很低、信号采样点数很少等情况下，相比于已有的频率估计方法，本发明方法能够在正弦信号频率附近得到更加明显的谱峰，从而得到更加精确的频率估计值。

附图说明

图1为本发明所提供的基于自适应滤波器幅频响应的高精度频率估计方法的流程图；

图2为仿真实验中不同的自适应滤波器中心频率对应的自适应滤波器幅频响应；

图3-1为仿真实验中当信噪比SNR＝10dB，信号采样点数N＝2000时，FFT频率估计方法的功率谱。

图3-2为仿真实验中当信噪比SNR＝10dB，信号采样点数N＝2000时，MUSIC频率估计方法的功率谱。

图3-3为仿真实验中当信噪比SNR＝10dB，信号采样点数N＝2000时，本发明方法在第一个正弦信号频率附近的伪功率谱。

图3-4为仿真实验中当信噪比SNR＝10dB，信号采样点数N＝2000时，本发明方法在第二个正弦信号频率附近的伪功率谱。

图4-1为仿真实验中当信噪比SNR＝-20dB，信号采样点数N＝500时，FFT频率估计方法的功率谱。

图4-2为仿真实验中当信噪比SNR＝-20dB，信号采样点数N＝500时，MUSIC频率估计方法的功率谱。

图4-3为仿真实验中当信噪比SNR＝-20dB，信号采样点数N＝500时，本发明方法在第一个正弦信号频率附近的伪功率谱。

图4-4为仿真实验中当信噪比SNR＝-20dB，信号采样点数N＝500时，本发明方法在第二个正弦信号频率附近的伪功率谱。

图5为仿真实验中在不同的信噪比下已有频率估计方法和本发明方法的频率估计均方根误差。

图6为仿真实验中在不同的信号采样点数下已有频率估计方法和本发明方法的频率估计均方根误差。

具体实施方式

下面结合附图对本发明提供的基于自适应滤波器幅频响应的高精度频率估计方法进行详细说明。

图1是本发明所提供的基于自适应滤波器幅频响应的高精度频率估计方法的流程图。该流程图的第一步是计算信号协方差矩阵，首先，对观测的连续时间信号进行采样得到离散时间信号，然后，由离散时间信号计算得到信号协方差矩阵R。第二步是确定第一个正弦信号的频率区间，即由先验信息确定第一个正弦信号的频率区间Φ₁。第三步是计算自适应滤波器的权系数向量，即由信号协方差矩阵R和频率区间Φ₁计算自适应滤波器的权系数向量h。第四步是计算自适应滤波器的幅频响应，即由自适应滤波器的权系数向量h计算自适应滤波器在频率区间Φ₁内的幅频响应H(f)。第五步是计算信号伪功率谱，即由自适应滤波器的幅频响应在频率区间内的最小值及最大值计算信号伪功率谱P(f_c)。第六步是计算第一个正弦信号的频率估计值，即在频率区间Φ₁内搜索信号伪功率谱P(f_c)的谱峰，并将信号伪功率谱P(f_c)的谱峰对应的频率作为第一个正弦信号的频率估计值f₁。若剩余未估计的正弦信号个数大于0，则根据先验信息确定下一个正弦信号的频率区间，然后重复执行步骤三至步骤六，得到该正弦信号的频率，以此类推，最终得到所有正弦信号的频率。

图2至图6是利用本发明方法在具体实施方式下进行仿真实验的处理结果，其中仿真实验的基本参数设置如下：

观测信号包含正弦信号和噪声，其中正弦信号的个数K＝2，正弦信号的频率分别为f₁＝150Hz和f₂＝300Hz，噪声为零均值加性高斯白噪声。对观测信号的采样频率为1000Hz，信号向量x(k)的长度m取为m＝100。在FFT方法中，设置FFT运算点数为65536点，在MUSIC方法、Root-MUSIC方法、ESPRIT方法以及本发明方法中，设置正弦信号的个数为2，在本发明方法中先验信息还包括：第一个正弦信号的频率区间为Φ₁＝[145Hz,155Hz]，第二个正弦信号的频率区间为Φ₂＝[295Hz,305Hz]。

图2是仿真实验中当信噪比SNR＝10dB，信号采样点数N＝2000时，不同的自适应滤波器中心频率对应的自适应滤波器幅频响应。图2中横轴表示频率，单位为Hz，纵轴表示自适应滤波器幅频响应的大小。在图2中，当自适应滤波器中心频率f_c＝150Hz时，即自适应滤波器中心频率f_c等于第一个正弦信号的频率f₁时，自适应滤波器幅频响应在频率区间Φ₁内近似为常数，且不存在明显的零陷。当自适应滤波器中心频率f_c＝149.9Hz或f_c＝149.8Hz时，即自适应滤波器中心频率f_c不等于第一个正弦信号的频率f₁，此时自适应滤波器幅频响应在频率区间Φ₁内存在明显的零陷。

图3-1至图3-4是仿真实验中当信噪比SNR＝10dB，信号采样点数N＝2000时，FFT频率估计方法和MUSIC频率估计方法的功率谱以及本发明方法在两个正弦信号频率附近的伪功率谱。在图3-1至图3-4中，横轴表示频率，单位为Hz，纵轴表示归一化功率谱或归一化伪功率谱，单位为dB。由图3-1至图3-4可见，FFT方法、MUSIC方法以及本发明方法在正弦信号的频率附近均能得到明显的谱峰，因此，FFT方法、MUSIC方法以及本发明方法在信噪比足够大和信号采样点数足够多的条件下，都能够得到精确的频率估计值。

图4-1至图4-4是仿真实验中当信噪比SNR＝-20dB，信号采样点数N＝500时，FFT频率估计方法和MUSIC频率估计方法的功率谱以及本发明方法在两个正弦信号频率附近的伪功率谱。在图4-1至图4-4中，横轴表示频率，单位为Hz，纵轴表示归一化功率谱或归一化伪功率谱，单位为dB。由图4-1至图4-4可见，由于信噪比很低且信号采样点数很少，由FFT方法计算得到的功率谱没有得到明显的谱峰，而由MUSIC方法计算得到的功率谱在错误的频率上得到了谱峰。但是，本发明方法仍然能够在正弦信号的频率附近得到较明显的谱峰。因此，相对于FFT方法和MUSIC方法等已有频率估计方法，本发明方法在信噪比很低和信号采样点数很少等情况下能够取得更加精确的频率估计。

图5为仿真实验中当信号采样点数N＝500时FFT方法、MUSIC方法、Root-MUSIC方法及ESPRIT方法等已有频率估计方法和本发明方法在不同信噪比下的频率估计均方根误差。图5中横轴表示信噪比，单位为dB，信噪比的变化范围是-30dB～10dB，纵轴表示频率估计均方根误差，单位为Hz。由图5可见，本发明方法的频率估计精度在不同信噪比下均优于或接近已有频率估计方法，特别是当信噪比很低时，例如当信噪比SNR＜-10dB时，本发明方法的频率估计均方根误差明显低于已有频率估计方法。

图6为仿真实验中当信噪比SNR＝-20dB时FFT方法、MUSIC方法、Root-MUSIC方法及ESPRIT方法等已有频率估计方法和本发明方法在不同的信号采样点数下的频率估计均方根误差。图6中横轴表示信号采样点数，信号采样点数的变化范围是100～1000，纵轴表示频率估计均方根误差，单位为Hz。由图6可见，本发明方法的频率估计精度在不同信号采样点数下均优于已有频率估计方法。

Claims (1)

1.一种基于自适应滤波器幅频响应的高精度频率估计方法，其特征在于，所述方法包括下述步骤：

第一步：计算信号协方差矩阵

首先，对观测的连续时间信号进行采样得到离散时间信号{x(1),x(2),...,x(N)}，其中采样频率为f_s，单位为Hz，N表示离散时间信号的采样点数，然后，由离散时间信号{x(1),x(2),...,x(N)}构造信号向量x(k)，即：

x(k)＝[x(k-m+1),x(k-m+2),...,x(k)]^H

其中，m为信号向量x(k)的长度，m为正整数，且1＜m＜N，k也为正整数，且m≤k≤N，m的具体大小根据实际应用情况确定，符号[·]^H表示共轭转置操作；最后，由信号向量x(k)计算信号协方差矩阵R，即：

$R=\frac{1}{N}\underset{k=m}{\overset{N}{\Σ}}x\left(k\right)x^H\left(k\right)$

其中x^H(k)表示对信号向量x(k)的共轭转置；

第二步：确定第一个正弦信号的频率区间

根据先验信息确定第一个正弦信号的频率区间，并记为Φ₁＝[f₁ ^l,f₁ ^r]，其中f₁ ^l为频率区间的下边界，f₁ ^r为频率区间的上边界，单位均为Hz；第一个正弦信号的频率落在频率区间Φ₁之中；

第三步：计算自适应滤波器的权系数向量

根据第一步得到的信号协方差矩阵R以及第二步得到的第一个正弦信号的频率区间Φ₁，按照下式计算自适应滤波器的权系数向量h，即：

$h=\frac{R^{-1}a\left(f_c\right)}{a^H\left(f_c\right)R^{-1}a\left(f_c\right)},f_c\∈{\mathrm{\Φ}}_1$

其中，R^-1表示信号协方差矩阵R的逆，向量f_c为自适应滤波器的中心频率，f_c∈Φ₁表示f_c的取值范围在频率区间Φ₁内；按照上述方式设计的自适应滤波器具有以下特性：当自适应滤波器的中心频率不等于第一个正弦信号的频率时，自适应滤波器的幅频响应在频率区间内存在明显的零陷，而当自适应滤波器的中心频率等于第一个正弦信号的频率时，自适应滤波器的幅频响应在频率区间内不存在明显的零陷；

第四步：计算自适应滤波器的幅频响应

根据第三步计算得到的自适应滤波器的权系数向量h，计算自适应滤波器在频率区间Φ₁内的幅频响应H(f)，即：

H(f)＝|a^H(f)h|,f∈Φ₁

其中，向量|·|表示求复数的幅度；自适应滤波器在频率区间Φ₁内的幅频响应H(f)也可以由信号协方差矩阵R和自适应滤波器的中心频率f_c直接计算得到，即：

$H\left(f\right)=\left|\frac{a^H\left(f\right)R^{-1}a\left(f_c\right)}{a^H\left(f_c\right)R^{-1}a\left(f_c\right)}\right|,f\∈{\mathrm{\Φ}}_1$

第五步：计算信号伪功率谱

首先，根据第四步计算得到的自适应滤波器在频率区间Φ₁内的幅频响应H(f)，搜索其最大值和最小值，分别记为：

$H_{max}\left(f_c\right)=\underset{f\∈{\mathrm{\Φ}}_1}{max}\left|\frac{a^H\left(f\right)R^{-1}a\left(f_c\right)}{a^H\left(f_c\right)R^{-1}a\left(f_c\right)}\right|$

$H_{\min }\left(f_c\right)=\underset{f\∈{\mathrm{\Φ}}_1}{min}\left|\frac{a^H\left(f\right)R^{-1}a\left(f_c\right)}{a^H\left(f_c\right)R^{-1}a\left(f_c\right)}\right|$

然后，由自适应滤波器幅频响应最大值H_max(f_c)和自适应滤波器幅频响应最小值H_min(f_c)计算信号伪功率谱P(f_c)，即：

$P\left(f_c\right)=H_{\min }^2\left(f_c\right)/H_{\max }^2\left(f_c\right)$

信号伪功率谱P(f_c)也可以由信号协方差矩阵R和自适应滤波器的中心频率f_c直接计算得到，即：

$P\left(f_c\right)={\left\{\underset{f\∈{\mathrm{\Φ}}_1}{\min }\left|a^H\left(f\right)R^{-1}a\left(f_c\right)\right|/\underset{f\∈{\mathrm{\Φ}}_1}{\max }\left|a^H\left(f\right)R^{-1}a\left(f_c\right)\right|\right\}}^2$

第六步：计算第一个正弦信号的频率估计值

根据第五步计算得到的信号伪功率谱P(f_c)，搜索其在频率区间Φ₁内的谱峰，并将信号伪功率谱P(f_c)的谱峰对应的频率作为第一个正弦信号的频率估计值f₁，即：

$f_1=\arg \underset{f_c\∈{\mathrm{\Φ}}_1}{max}{\left\{\underset{f\∈{\mathrm{\Φ}}_1}{\min }\right|a^H\left(f\right)R^{-1}a\left(f_c\right)|/\underset{f\∈{\mathrm{\Φ}}_1}{max}|a^H\left(f\right)R^{-1}a\left(f_c\right)\left|\right\}}^2$

若在观测信号中正弦信号的个数大于1，则首先根据先验信息确定其它正弦信号的频率区间，然后重复上述步骤中的第三步至第六步，即可得到其它正弦信号的频率估计值。