项目(9)-音频信号分析仪(快速傅里叶变换)

音频信号分析仪

一、项目简介

信号频谱分析是信号与系统和数字信号处理里面的重要内容。它在医学、通信、语音、图像等众多领域上有着重大作用和意义。音频信号是我们日常生活接触最多的信号之一，许多的音频处理都需要建立在音频频谱分析的基础上。该项目由增强型STC12LE5A60S2单片机为主控制器，通过单片机内部AD转换，对音频信号进行采样，把连续信号离散化，然后通过FFT快速傅里叶变换运算，在时域和频域对音频信号各个频率分量以及功率等指标进行分析和处理，最后通过12864液晶对信号的频谱进行显示。

二、项目要求

1、输入信号频率范围20Hz-10KHz。

2、FFT的分辨能力尽量高(即FFT点数尽量大)。

3、能得出信号最大频率分量的频率。

4、用单片机产生多个自测信号供测试用。

5、扩展要求1：用Rife算法进行频率校正。

6、扩展要求2：在液晶上显示信号的动态频谱图。

7、扩展要求3：输入信号频率范围20Hz-20KHz。

三、项目方案

AD转换：开发板选用的处理器是STC12LE5A60S2，其内部自带10位8通道AD转换器，因此不需要自行设计AD转换电路，直接使用内部AD对输入信号采样即可。

FFT点数：单片机内部程序存储空间(flash)为60K字节，而数据存储空间(SRAM)仅有1280字节。由于FFT浮点运算对单片机的RAM和速度的要求都比较高，而且按FFT运算的特点，FFT点数越大，分辨率越高，消耗的RAM空间也越大。假设FFT点数为N，那么可以计算出FFT运算实部XR[N]和虚部XR[N]所消耗的空间为2*4*N字节，即8N。而基2FFT蝶形运算时N必须是2的指数次。由此可知，N最大只能取128。这时实部和虚部共占用1024字节。但是，整个程序除了FFT运算，还有很多地方要消耗不少的RAM空间，比如液晶显示部分，各种函数或变量等等。因此如果取N=128的话，单片机的RAM空间就会变得很紧张，这对项目开发非常不利。为了保险起见，先取FFT的点数N=64。

采样率与分辨率：根据奈奎斯特定理，采样率>=两倍信号最高频率。基本要求中信号频率范围是20Hz-10KHz，而扩展要求的信号频率范围是20Hz-20KHz。为了达到扩展要求，应该将采样率设置成40KHz或以上。分辨率=采样率fs/点数N，点数N已确定是64，那么采样率越小，分辨率的值越小，也就是能分辨的最小频率间隔刻度越小，这就越精确。因此该项目选择40KHz采样率。这样在保证满足奈奎斯特定理的条件下能达到最高分辨率。采样时刻基准由定时器溢出中断提供。

自测信号：可用单片机内部PWM产生方波信号，也可用定时器来产生方波信号。满足信号频率范围20Hz-20KHz即可。

信号经采样和FFT变换后，得出信号的频谱或功率谱，直接输出到液晶上显示，即可得到信号的频谱图。可用Rife算法对FFT求得的最大频率值进行校正，最后输出显示。

四、相关原理知识

4.1、FFT原理简介

FFT，即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。设x(n)为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）, 即N点DFT变换大约就需要N^2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列（设N=2k,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）^2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2（N/2）^2=N+N^2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N 点的DFT变换就只需要Nlog(2)(N)次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT 的优越性。

4.2、离散傅里叶变换DFT

为确定输入信号采样的频谱，我们需要对这些输入采样进行离散傅立叶变换(DFT)。DFT的定义如下：

其中N是采样的数量，X(k)是频谱，x(n)是一组输入采样。利用欧拉等式展开求和符，并分离输入采样和频谱的实部和虚部，得到以下等式：

式2和3中，求和符中第二项的消失是由于输入采样全部为实数。假定我们有N个采样，直接计算式2和3需要2N²次乘法和2N(N - 1)次加法。这样，我们的256点输入采样DFT将需要进行131,072次乘法和130,560次加法运算。

4.3、radix-2 FFT算法

有多种FFT算法可供使用。该项采用普通的radix-2算法，继续将DFT分解为两个更小的DFT。为此，N必须是2的指数。这种radix-2 FFT算法的步骤可归纳如图2所示的蝶型运算。观察这些蝶型运算我们可以发现，radix-2算法仅需(N / 2)log2(N)次乘法和Nlog2(N)次加法。图2中用到的参数WN就是通常所谓的“旋转因子”，可以在执行算法前预先计算出来。

图2. 利用蝶型运算实现N = 8的FFT。

在图2中，FFT的输入显示为一种特殊的排列顺序，这种序列是对原始序列索引号的二进制位反转后得到的。因此，当我们对N = 8个采样执行radix-2 FFT 算法时，需要将输入数据的原始序列：

0 (000b), 1 (001b), 2 (010b), 3 (011b), 4 (100b), 5(101b), 6(110b), 7(111b)

重新排列为：

0 (000b), 4, (100b), 2 (010b), 6 (110b), 1 (001b), 5 (101), 3 (011), 7 (111)

FFT输出则以正确的顺序排列。图2还说明，每个单独的蝶型运算所得的结果，是下一级FFT运算所需的唯一数据。由于运算过程可“即位”进行，新值可替代旧值，这样，计算N个采样的FFT只需要2N个变量(因为每个数据都包括实部和虚部两部分)。

FFT完成后，结果为复数形式。式4和5将结果转换为极坐标方式后表示为：

有关DSP的文献中可以找到很多优化方法，可使上述DFT/FFT算法更小或更快。其中最重要的一种优化方法(可能也是最容易实现的)源于这样一个事实，那就是作为一个实数信号，其DFT幅度是相关于X(N / 2)对称的，因此：

4.4、Rife算法频率校正

若是{g(n)}经FFT变换后功率谱的最大值，Rife算法给出的频率估计公式为

式中：，当时，，否则为-1；fs为采样频率；k0为

功率谱最大值对应的量化频率，k0为整数。修正因子0<=<=1/2。根据功率谱的对称性特点，被估计频率f0介于k0fs/N与(k0+1/2α)fs/N之间。

f0在最大频率谱线与次大频率谱线之间，利用rife算法可把正弦频率校正成较准确。也可利用插值算法进行频率校正，但由于该算法太过繁复，运算量较大，因此不采用。

4.6、频谱泄露与加窗处理