FFT三点插值法估计频率

zwchhm

已于 2024-02-22 17:05:29 修改

阅读量1k

点赞数 9

文章标签：算法

于 2024-02-22 17:04:19 首次发布

本文链接：https://blog.csdn.net/zwchhm/article/details/136236649

版权

本文介绍了在FFT频率估计中，通过三点插值法提高精度的方法，包括插值思路、公式推导，以及使用1024点正弦波的仿真示例。核心内容是利用抛物线模型求解插值参数，以补偿FFT精度不足的问题。

摘要由CSDN通过智能技术生成

前言

利用FFT进行频率的估计时精度较低，若信号采样率为 $f_s$ ，FFT点数为 $N_s$ ，则信号精度仅为 $\delta_F = f_s / N_s$ ，需要对FFT后的频率数据进行插值来获取更高的精度。本文推导三点插值公式。

一、插值思路

可用FFT后的信号模值或者模值的对数进行插值。以信号模值为例，设模值最大点为 $y_m$ ，对应的频点为 $x_m$ ，左右两边点的模值分别为 $y_l$ 和 $y_r$ ，对应的频点分别为 $x_l$ 和 $x_r$ 。则三点插值的思路是将（对数）模值波形近似为抛物线，利用抛物线公式找出最大值。

设抛物线公式为

$y = ax^2 + bx + c$

插值思路就是利用三点求出抛物线系数a、b和c的值，并计算得到极值点对应的频点

$x^* = -\frac{b}{2a}$

二、插值公式推导

将三点代入公式可得：

$y_l = ax_l^2 + bx_l + c$ （1）

$y_m = ax_m^2 + bx_m + c$ (2)

$y_r = ax_r^2 + bx_r + c$ (3)

且满足 $x_m - x_l = x_r - x_m = \delta_f$ ，即 $x_l$ 、 $x_m$ 和 $x_r$ 是一个等差数列。

由公式（3）-公式（1）可得：

$y_r - y_l = a(x_r + x_l)(x_r - x_l) + b(x_r - x_l)$

整理可得：

$b = \frac{y_r - y_l}{2\delta_f} - 2ax_m$ (4)

再由公式（1）+公式（3）-2公式（2）可得：

$y_l + y_r - 2y_m = a(x_l^2 - x_m^2 + x_r^2 - x_m^2) + b(x_l + x_r - 2x_m)$ (5)

根据 $x_l$ 、 $x_m$ 和 $x_r$ 等差数列的关系可得：

$x_l^2 - x_m^2 + x_r^2 - x_m^2 = -\delta_f(x_l + x_m) + \delta_f(x_r + x_m) = 2\delta_f^2$

$x_l + x_r - 2x_m = 0$

代入公式（5）可解出a为

$a = \frac{y_l + y_r - 2y_m}{2\delta_f^2}$ (6)

由公式（4）和公式（6）可得极值点对应的频偏为

$x^* = -\frac{b}{2a} = -\frac{\frac{y_r - y_l}{2\delta_f} - 2ax_m}{2a} = x_m + \frac{2(y_l - y_r)}{y_l + y_r - 2y_m}\delta_f$

上式中，第一项 $x_m$ 是FFT后最大值点对应的频点，第二项是通过插值需要补偿的频偏。

三、仿真验证

以1024点正弦波为例，进行1024点的FFT，插值核心代码如下：

function fe = FFT_Interpolate(FFT_in, nfft)
% FFT_in -- 复数正弦波
% nfft -- FFT点数
Freq_Domain = fft(FFT_in, nfft);
abs_F = abs(Freq_Domain);
[max_value, index] = max(abs_F);
if(index == 1)
    left_value = abs_F(nfft);
    right_value = abs_F(2);
elseif(index == nfft)
    left_value = abs_F(nfft - 1);
    right_value = abs_F(1);
else
    left_value = abs_F(index - 1);
    right_value = abs_F(index + 1);
end


delta_f = (left_value - right_value) / (right_value + left_value - 2 * max_value) / 2;
fe = (index - 1 + delta_f) * fs / nfft;
if(fe > fs / 2)
    fe = fe - fs;
end