马旭峰，邢阳阳，周子程，范雯，尹鸿峰
（沧州交通学院 计算机与信息技术学院，河北黄骅，061100）
0 引言
信号处理技术已经深入到各个领域中去，对信号的分析也早已从经典的时域分析方法转变为频域分析方法，对于音频信号同样如此。音频信号主要分为两类：语音信号和非语音信号。非语音信号具有非常简单的语意和语法信息，信息量低且能简单地识别，而语音信号则更加复杂多变，对其进行时域分析往往不能发现更深层次的信息，但是它的频谱对外部环境的改变则表现出一定的顽健性，有很明显的声学特性。
频域分析的核心理论为离散傅里叶变换（Discrete Fourier Transform，DFT），是信号分析最常用的工具[1]，然而DFT只是对信号进行理论上的分析，无法运用到实践中，最主要的原因在于这种方法的计算量非常大。本文采用了它的快速实现算法，即快速傅里叶变换（FFT），此算法充分地利用了DFT运算中的两个主要性质：对称性和周期性。极大地减小了运算量，计算效率提高了近100倍。
采用DSP微处理器作为该系统的核心单元，利用其在数字信号处理领域高性能、低功耗优点，结合一款语音编解码芯片TLV320AIC23B（以下简称AIC23B）完成系统整体设计。
本文主要贡献在于设计了一种基于DSP的频谱分析系统，可以实现对音频信号的非实时的频谱分析。
1 系统框架
频谱分析系统的硬件实现方案如图 1所示。其中，TMS320VC5502型号定点DSP处理器为系统的核心，它的部分技术参数如下[3～6]：
图1 系统硬件框架
(1)300MHz处理器主频；
(2)2M*32位SDRAM，工作时钟频率为100MHz；
(3)16K*8位的EEPROM，IIC串行接口，串行位时钟400Kbps；
(4)256K*16位Flash。
音频信号的采集与转换使用的是型号为AIC23B的语音芯片，该芯片内部采用了Sigma-delta 过采样技术，可以在8K～96K采样率范围内提供多种位数选择的抽样技术。AIC23B控制接口具有I2C和SPI两种工作方式，由于DSP芯片内具有I2C总线功能，因此配置AIC23B控制接口为I2C模式，实现了DSP处理器对语音芯片的控制通道。而AIC23B的数据接口主要的工作方式也是两种，即IIS模式和DSP模式。本次设计配置AIC23B数据接口为DSP模式，利用DSP芯片内部McBSP串口连接AIC23B的数据接口进行数据传输。
2 频谱分析算法原理
离散傅里叶变换简称DFT（Discrete Fourier Trans form）在数字信号处理技术中占据着极为广泛的地位，它是处理有限长序列的有效方法。
当X(n)为一个有限长序列时，其离散傅里叶正反变换[2]为：
对于x(n)和X(k)，知道二者任意一个值，就能计算出另一个。因此，可以给它们乘以相应的内插函数来确定连续函数，因而离散傅里叶变换可以看作是连续傅里叶的近似。
但在很长一段时间内，我们无法将DFT实现大范围的推广，这是因为DFT计算起来相当繁琐，无法高效地应用到计算机中。直到快速傅里叶变换（FFT）算法的推出和计算机技术的进步，极大程度上压缩了信号处理的计算时间，从而使得DFT在许多工程领域中得到普遍的使用。本文选择的是按频率抽选的基-2FFT方法。
设序列的点数为N=2L，L为整数。先将输入部分按N的序号分成上下两个子序列。
按照k的奇偶顺序将X(k)分为两部分。式(4)可变为：
令：
则式(5)可变为：
其中，上式(6)所反映的计算关系可由图 2中的蝶形来形象示出。
图2 按频率抽选的蝶形运算示意图
由图可知，每一次蝶形运算都需要进行一次乘法和两次加法运算。这样就可以将一个N点DFT按k的奇偶分成了2个N/2点的DFT了。假如当N=8时，上述的分解过程如图 3所示。
图3 一个N=8点DFT分解为2个N/2点DFT运算示意图
由于分解完成的N/2的值还是偶数，因此可继续将N/2点的DFT的结果按照奇偶分开为两部分，这就将N/2点的DFT继续分成为2个N/4点的DFT。这样的分解得到的最终结果就只有加减操作，极大程度地减少了计算量。
3 频谱分析算法实现与改进
3.1 算法实现
本次用到的频谱分析方法是基2按频率抽选的快速傅里叶变换。FFT算法程序通过三层循环完成[1]。
第二层循环：完成系数r NW的运算，也就是求旋转因子中r。
本次使用的DSP芯片的型号属于16位定点DSP，而本文要实现的运算属于浮点型的运算。为了实现运算结果并达到计算精度，需要对运算的数据进行定标操作，也就是在编写程序的时候应该对小数点的位置进行确定。本次使用的定标方法是Q格式定标法。
具体做法是：先对所有旋转因子都放大2^Q倍，为保证旋转因子的差异化，Q必须大于等于L，旋转因子被放大后，为了保证其模为1，要在最内层循环中把每一次蝶形运算的结果右移Q位来抵消旋转因子放大，从而得到正确的结果。通过简单的右移位运算代替除法运算来抵消之前的放大的办法大大地节省了运算时间。
需要注意是Q的值越大，数据精确度变高，对器件要求会高。因此，对于定点数值，它的数据区间和精度是有冲突的。在实际的算法中，要充分考虑到Q值的选取，以达到最好的状态。本设计令Q的值为15，即旋转因子被放大了32767倍。
3.2 算法改进
本文对FFT中倒位序算法程序进行了初步改进。方法为使用了一种特殊点数的程序来代替通用的倒位序程序，具体代码如下：
这段倒位序程序代码是针对点数为512点的FFT变换，主要是通过移位运算来求该数的倒位序，相比于通用的倒位序程序，此代码效率更高，测试结果如下：
以512点为例，统计了FFT代码的效率，算法改进之前的统计结果如图 4所示。其中，INCL TOTAL指这段代码消耗得总的时钟周期数。
图4 改进之前FFT算法效率
由图4可以看出，512点FFT所消耗的时钟周期数为626720，与此同时main函数消耗的时钟周期数为3999184。而本次设计的CPU主频为240MHz（最高值为300MHz），可以计算得到FFT算法的代码运行时间为2.611ms。
而经过改进之后，如图5所示，FFT算法程序消耗的时钟周期数由原来的626720降低到601388，代码运行时间由2.611ms 减少到2.505ms，算法效率提高了4.23%。
图5 改进之后FFT算法效率
4 系统运行结果
将设备和仿真器连接，通过开发板上的Line插孔向语音芯片内输入音频，在CSS集成开发环境下可以看到变换结果。信号采集后的时域、幅值、功率、对数等图形结果如图6～图10所示。
图6 语音信号的时域图
图7 FFT变换图谱
图8 幅值谱
图9 功率谱
通过以上这几幅图形可以对该信号进行分析。其中，幅值图表示了该语音信号各个谐波分量的幅值随频率的线性分布情况。功率谱表示的是功率，也是语音信号各个谐波频率幅度的自乘，其反映出的是频率的成分。当要观察所有频率成分时，一般采用更接近于人耳实际听音时感觉的对数谱。
图10 对数谱
5 结论
针对各个工程领域对信号频谱分析的需求，设计了一种基于DSP的频谱分析系统，系统依靠了DSP微处理器在信号处理过程中运算高效的特点，采用C语言编写了各个寄存器配置程序以及FFT算法程序，改进并提高了FFT算法的执行效率。本系统目前只针对音频信号实现了对非实时采集的音频信号的频域变换，变换结果清晰地展示了信号的频域特征，具有直观的可视化效果。
