CN105307095B - 一种基于fft的高分辨率音频频率测量方法 - Google Patents

Abstract

本发明涉及音频频率测量技术领域，具体涉及一种基于FFT的高分辨率音频频率测量方法，包括，获取波形数据；对所述波形数据进行加汉宁窗快速傅里叶变换FFT；选取频谱中目标信号处的峰值；计算峰值左侧点FFT的结果和右侧点FFT的结果，并比较该左侧结果和右侧结果的大小；计算得到参数k，y₁，y₂；根据该参数k，y₁，y₂，计算所述波行数据的频率及幅度。本发明可以明显提高频率及幅度的计算结果精度，同时可以保证在分辨率精度较低时准确的识别出相邻频谱的频率及幅度；对于50kHz采样率、1024点FFT，利用该算法，频率分辨率可由常规算法的48Hz提升到0.1Hz。

Description

一种基于FFT的高分辨率音频频率测量方法

技术领域

本发明涉及音频测试技术领域，具体涉及一种基于FFT的高分辨率音频频率测量方法。

背景技术

音频是多媒体中的一种重要媒体，蕴含了丰富的语义，有效地丰富和补充了人们对信息的感知和获取。声音经过模拟设备记录或再生，成为模拟音频，再经数字化成为数字音频，音频分析就是以数字音频信号为基础，以数字信号处理为分析手段，提取音频信号在时域、频域内的一系列特性的过程。音频分析被广泛应用于测量各类音频系统的时域特性、频域特性及失真特性等。

在音频分析过程中，最常用的方法就是傅里叶变换法。快速傅里叶变换FFT是音频分析中的常用算法，快速傅里叶变换存在栅栏效应和频谱泄露现象，使算出的信号参数即频率、幅值和相位不够准确。对于传统的利用FFT进行音频分析的方法中，频率误差大小取决于实际频率值是否刚好落在N点快速傅里叶变换的对应频率通道上，如果实际频率落在两个相邻频率通道中间，那对应的误差就会比较大，且误差量不可控。同时，频率误差有大小还取决于FFT点数，FFT点数越小，直接计算出的频率误差也会进一步恶化，无法满足准确的频率测量要求。

FFT作为音频分析中的常用算法，是离散傅里叶变换(DFT)的一种快速算法，它通过将长序列的DFT分解为短序列的DFT进行计算，从而使运算量大大减少，利用快速傅里叶变换可计算得到有限长信号的离散频谱。

发明内容

本发明的目的是提供一种基于FFT的高分辨率音频频率测量方法，提高频率及幅度的计算结果精度。

为达上述目的，本发明提供了一种基于FFT的高分辨率音频频率测量方法，包括：

获取波形数据；

对所述波形数据进行加汉宁窗快速傅里叶变换FFT；

选取频谱中目标信号处的峰值；

计算峰值左侧点FFT的结果和右侧点FFT的结果，并比较该左侧结果和右侧结果的大小；

计算得到参数k，y₁，y₂；

根据该参数k，y₁，y₂，计算所述波行数据的频率及幅度。

进一步的，

汉宁窗函数为：

W(n)＝0.5-0.5*cos(2π*n/N)n＝0～N-1；

所述对所述波形数据进行加汉宁窗快速傅里叶变换FFT具体为：

fftbuf[n]＝FFT(wavebuf[n]*W(n))n＝0～N-1；

上式中，FFT(...)代表快速傅里叶变换运算过程，fftbuf为用于存放FFT运算结果的数组，wavebuf为采集到的原始波形数据，N为采样点数。

进一步的，

所述计算得到参数k，y₁，y₂，具体包括：

若fftbuf[maxI-1]>fftbuf[maxI+1]，则选取峰值及其左边的点进行计算，此时可得出以下参数：

k＝maxI-1

y₁＝|fftbuf[maxI-1]| y₂＝|fftbuf[maxI]|

其中，maxI表示所述目标信号处的峰值的位置；k代表峰值左边临近点的位置，y₁、y₂分别代表峰值左边临近点处及峰值处的FFT运算结果；

若fftbuf[maxI-1]≤fftbuf[maxI+1]，则选取最大值及其右边的点进行计算，此时可得出以下参数：

k＝maxI

y₁＝|fftbuf[maxI]| y₂＝|fftbuf[maxI+1]|

其中，k代表峰值右边临近点的位置，y₁、y₂分别代表峰值处及峰值右边临近点处的FFT运算结果。

进一步的，

所述根据该参数k，y₁，y₂计算所述波行数据的频率及幅度具备包括：

其中，

k₀＝α+k+0.5；

式中，freq、ampl分别为计算得到的频率结果以及幅度结果，fullfreq为采样率，N为采样点数，k₁、k₂、k₃、k₄为公式推导得出的具体系数。

本发明能够达到以下有益效果：

本发明在现有FFT运算方法的基础上，提出了一种改进的计算频谱频率及幅度的方法，通过：获取波形数据；对所述波形数据进行加汉宁窗快速傅里叶变换FFT；选取频谱中目标信号处的峰值；计算峰值左侧点FFT的结果和右侧点FFT的结果，并比较该左侧结果和右侧结果的大小；计算得到参数k，y₁，y₂；根据该参数k，y₁，y₂，计算所述波行数据的频率及幅度；该方法可以明显提高频率及幅度的计算结果精度，同时可以保证在分辨率精度较低时准确的识别出相邻频谱的频率及幅度；对于50kHz采样率、1024点FFT，利用该算法，频率分辨率可由常规算法的48Hz提升到0.1Hz。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明一种基于FFT的高分辨率音频频率测量方法的流程图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在实际进行数字信号处理时，都需要把信号的观察时间限制在一定的时间间隔内，只需要选择一段时间信号对其进行分析。也就是取用有限个数据，将信号截断的过程相当于将信号进行加窗函数的操作。经过这样的操作后，往往会发生频谱分量从其正常频谱扩展开来的现象，即称为“频谱泄露”现象。当利用传统的傅里叶变换进行音频分析时，时域中的截断是必需的，所以泄露效应也是快速傅里叶变换所存在的，这就会大大降低计算所得信号的频率、幅度等参数，必须对其进行抑制。

考虑到传统傅里叶变换存在栅栏效应和频谱泄露现象，使计算得到的频谱结果精度较低，现有的方案通常是先选择特定的窗函数对傅里叶变换进行加窗处理，设x(n)是一个长序列，w(n)是长度为N的窗函数，用w(n)截断x(n)，得到N点序列x_N(n)，形式如下：

x_N(n)＝x(n)*w(n)

在频域上则有，

然后，对加窗处理后的短序列进行快速傅里叶变换，可以得到离散频谱，即

式中，X(k)为频域信号，N为采样点数。

最后，在快速傅里叶变换得到的频谱中找到峰值所在位置，进而计算得到我们所需要的频率。

式中，freq为计算得到的频率，maxI为峰值所在位置，fullfreq为采样率，N为采样点数，f为频率分辨率。

如采样率为50KHz，采样点数为1024的FFT算法中，计算得到峰值所在位置为100，则计算可得频率值为4882Hz，频率分辨率为48.82Hz，可以看出频率分辨率较低。

当窗函数应用于FFT算法时，可以有效减少频谱泄露所引起的误差，选择不同的窗函数会得到不同的计算精度，为了尽可能的提高计算精度，首先就要选择一个频谱特性较好的窗函数，常用的窗函数有矩形窗、三角窗、海明窗、汉宁窗、布莱克曼窗、布莱克曼哈里斯窗，其中，加汉宁窗截断的正弦信号的频谱泄露主瓣的宽度为4个谱线间隔，要分辨出相邻的谐波，加汉宁窗差值FFT算法至少需要2个周期的采样点，所需要的采样周期数较少，可以很好的满足我们对高实时性及高计算精度的要求，所以在本发明中我们选取了汉宁窗对FFT运算进行加窗处理。

此外，当进行傅里叶变换时采样点数如果不够多，会导致计算得到的频谱峰值所在位置恰好处于两个采样点中间的情况，这会导致快速傅里叶变换得到的结果与实际值存在较大的偏差，利用传统的傅里叶变换得到的频谱在频率分辨率较低时，对相邻频谱的频率和幅度识别也不够精确，所以在本发明中，我们通过使用信号峰值处及其左右两侧的两个最临近点进行加权运算得到一个修正系数对频率进行修正以得到更精确的测量频率。

下面通过具体实例对本发明进行进一步阐释：

图1是本发明一种基于FFT的高分辨率音频频率测量方法的流程图。如图所示，该方法包括：

步骤101，获取波形数据；

步骤102，对所述波形数据进行加汉宁窗快速傅里叶变换FFT；

步骤103，选取频谱中目标信号处的峰值；

步骤104，计算峰值左侧点FFT的结果和右侧点FFT的结果，并比较该左侧结果和右侧结果的大小；

步骤105，计算得到参数k，y₁，y₂；

步骤106，根据该参数k，y₁，y₂，计算所述波行数据的频率及幅度。

进一步的，

汉宁窗函数为：

W(n)＝0.5-0.5*cos(2π*n/N)n＝0～N-1；

所述对所述波形数据进行加汉宁窗快速傅里叶变换FFT具体为：

fftbuf[n]＝FFT(wavebuf[n]*W(n))n＝0～N-1；

上式中，FFT(...)代表快速傅里叶变换运算过程，fftbuf为用于存放FFT运算结果的数组，wavebuf为采集到的原始波形数据，N为采样点数。

进一步的，

所述计算得到参数k，y₁，y₂，具体包括：

若fftbuf[maxI-1]>fftbuf[maxI+1]，则选取峰值及其左边的点进行计算，此时可得出以下参数：

k＝maxI-1

y₁＝|fftbuf[maxI-1]| y₂＝|fftbuf[maxI]|

其中，maxI表示所述目标信号处的峰值的位置；k代表峰值左边临近点的位置，y₁、y₂分别代表峰值左边临近点处及峰值处的FFT运算结果；

若fftbuf[maxI-1]≤fftbuf[maxI+1]，则选取最大值及其右边的点进行计算，此时可得出以下参数：

k＝maxI

y₁＝|fftbuf[maxI]| y₂＝|fftbuf[maxI+1]|

其中，k代表峰值右边临近点的位置，y₁、y₂分别代表峰值处及峰值右边临近点处的FFT运算结果。

进一步的，

所述根据该参数k，y₁，y₂计算所述波行数据的频率及幅度具备包括：

其中，

k₀＝α+k+0.5；

式中，freq、ampl分别为计算得到的频率结果以及幅度结果，fullfreq为采样率，N为采样点数，k₁、k₂、k₃、k₄为公式推导得出的具体系数。

本发明是一种基于FFT的高分辨率音频频率测量方法，其关键在于提高频率分辨率，现在常用方法是经过FFT计算便可得到信号的频率与幅度，经过本改进方法处理以后，利用传统FFT算法找到目标信号的峰值，然后选取该峰值左右两侧临近点的运算结果进行加权运算，可明显提高频率及幅度的计算结果精度；本发明的关键点之一在于用传统FFT算法找到目标信号的峰值后，选取该峰值左右两个最临近点中较大者与该峰值进行加权运算，具体的较大值选取方法可参照上文所述求取k、y1、y2的过程，加权运算的过程可参照计算参数及k0的过程；本发明应用于音频频率测量之后，可有效提高频率及幅度的计算结果精度，改进后的频率及幅度计算结果可由上文描述的、计算过程计算得到，利用该改进方法得到的频率及幅度相比于传统算法得到的频率及频率分辨率有了明显的提高。

通过以上实施例，本发明能够达到以下有益效果：

本技术在现有FFT运算方法的基础上，提出了一种改进的计算频谱频率及幅度的方法，本改进方法可以明显提高频率及幅度的计算结果精度，同时可以保证在分辨率精度较低时准确的识别出相邻频谱的频率及幅度；对于50kHz采样率、1024点FFT，利用该算法，频率分辨率可由常规算法的48Hz提升到0.1Hz。

本领域技术人员还可以了解到本发明实施例列出的各种说明性逻辑块(illustrative logical block)，单元，和步骤可以通过电子硬件、电脑软件，或两者的结合进行实现。为清楚展示硬件和软件的可替换性(interchangeability)，上述的各种说明性部件(illustrative components)，单元和步骤已经通用地描述了它们的功能。这样的功能是通过硬件还是软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护的范围。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种基于FFT的高分辨率音频频率测量方法，其特征在于，包括：

获取波形数据；

对所述波形数据进行加汉宁窗快速傅里叶变换FFT；

选取频谱中目标信号处的峰值；

计算峰值左侧点FFT的结果和右侧点FFT的结果，并比较该左侧结果和右侧结果的大小；

计算得到参数k，y₁，y₂；

根据该参数k，y₁，y₂，计算所述波形数据的频率及幅度；

所述计算得到参数k，y₁，y₂，具体包括：

若fftbuf[maxI-1]>fftbuf[maxI+1]，则选取峰值及其左边的点进行计算，此时可得出以下参数：

k＝maxI-1

y₁＝|fftbuf[maxI-1]| y₂＝|fftbuf[maxI]|

其中，maxI表示所述目标信号处的峰值的位置；k代表峰值左边临近点的位置，y₁，y₂分别代表峰值左边临近点处及峰值处的FFT运算结果，fftbuf为用于存放FFT运算结果的数组；

若fftbuf[maxI-1]≤fftbuf[maxI+1]，则选取最大值及其右边的点进行计算，此时可得出以下参数：

k＝maxI

y₁＝|fftbuf[maxI]| y₂＝|fftbuf[maxI+1]|

其中，k代表所述目标信号处的峰值的位置，y₁，y₂分别代表峰值处及峰值右边临近点处的FFT运算结果；

所述根据该参数k，y₁，y₂计算所述波形数据的频率：

fftbuf[maxI-1]>fftbuf[maxI+1]时，k₀＝α+k+0.5；

fftbuf[maxI-1]≤fftbuf[maxI+1]时，k₀＝α+k+0.5；

式中，freq为计算得到的频率结果，fullfreq为采样率，N为采样点数。

2.根据权利要求1所述的方法，其特征在于：

汉宁窗函数为：

W(n)＝0.5-0.5*cos(2π*n/N)，n＝0～N-1；

所述对所述波形数据进行加汉宁窗快速傅里叶变换FFT具体为：

fftbuf[n]＝FFT(wavebuf[n]*W(n))，n＝0～N-1；

上式中，FFT(...)代表快速傅里叶变换运算过程，fftbuf为用于存放FFT运算结果的数组，wavebuf为采集到的原始波形数据，N为采样点数。