STM32F103 如何实现 FFT?

[导读]笔者能力有限，如果文中出现错误的地方，还希望各位朋友能够给我指出来，我将不胜感激，谢谢~ 引言数字信号在我们生活中随处可见，自然而然地就会涉及到对于数字信号的处理，最为典型的一个应用就是示波器，在使用示波器的过程当中，我们会通过示波器测量到

笔者能力有限，如果文中出现错误的地方，还希望各位朋友能够给我指出来，我将不胜感激，谢谢~

引言

数字信号在我们生活中随处可见，自然而然地就会涉及到对于数字信号的处理，最为典型的一个应用就是示波器，在使用示波器的过程当中，我们会通过示波器测量到信号的频率以及幅值，同时我们也可以通过示波器对测量到的信号进行 FFT ，从而能够观察到待测信号的频谱，方便直观的看出信号的高频分量和低频分量，从而帮助我们去除信号中携带的噪声。而在嵌入式方面的应用，我们可以直接使用 DSP 芯片对信号进行处理，同时， ARM 公司推出的 Cortex-M4F 内核是带有 FPU ,DSP 和 SIMD 单元的，针对于这些单元也增加了专用的指令，指令如下图所示：

不同架构的指令集合

ARM 官方也对此专门做了一个 DSP 方面的库，方便用户调用。关于 Cortex M4 的信号处理本文暂不进行阐述，相反本文的对象是 Cortex M3 ,基于 STM32F103 的 FFT，而在上述图中，我们看到针对于 Cortex M3 来说，是不带 FPU 以及 DSP 的，那有如何来进行 FFT 呢？

FFT 的提出

在数字信号处理中常常需要使用到离散傅里叶变换(DFT)，从而能够获取到信号的频域特征。尽管传统的 DFT 算法能够获取到信号的频域特征，但是算法计算量大，耗时长，不利于进行计算机实时对信号进行处理。因此才有了 FFT 的出现。需要强调的是，FFT 并不是一种新的频域特征获取方式，而是 DFT 的一种快速实现算法。FFT 之所以能够改善运算量，是因为其充分运用了 DFT 运算中的对称性和周期性，从而能够将运算量从 N^2 减少到 N*log2(N),其中 N 为待计算的序列的长度。当 N 非常大的时候，这种优化在时间维度上的提升是非常显著的。下面是关于 DFT 和 FFT 所需乘法次数的比较曲线。

FFT 算法与 DFT 算法的比较

FFT 变换之后和原始信号的对应关系

假设我们对一个波形进行了采样，采样了 N 个点，经过 FFT 之后，就可以得到 N 个点的 FFT 结果，每一个点就对应着一个频率点。这个点的模值，就是该频率下的幅度特性。具体的关系就是如果原始信号的峰值为 A ，那么 FFT 的结果的每个点的模值就是 A 的 N / 2 倍。而第一个点就是直流分量，它的模值是直流分量的 N 倍。而每个点的相位就是在该频率下的信号的相位，第一个点表示的是直流分量，也就是 0 HZ的点，而最后一个点 N 的再下一个点（实际这个点是不存在的），也就是 N+1 个点则表示的是采样频率 Fs，这中间被 N - 1 个点平均分成 N 等份，每一个点的频率依次增加。也就是如果要计算某个点的频率，那么就只需要这样计算即可：Fn = (n - 1) * Fs / N。
从上述所展示的公式，我们可以知道 Fn 所能够分辨的频率为 Fs / N，如果采样频率 Fs 为 1024Hz,采样点数为 1024 点，则可以分辨到 1 HZ。也就是说采样 1s 时间的信号并做 FFT ，则结果可以分析精确到 1 Hz，如果采样 2 s 时间的信号并做 FFT，则结果可以分析精确到 0.5 Hz，所以也就说明了一个道理，如果要提高频率分辨率，则必须增加采样点数，也就是采样时间，下面这张图更能够清晰地表示这种关系：

频率分辨率

将原信号变换之后的频谱的宽度与原始信号也存在一定的关系。根据 Nyquist采样定律，FFT 之后的频谱宽度最大只能是原始信号采样率的 1/2，如果原始信号的采样频率为 4GS/s，那么 FFT 之后的频宽最多只能是 2GHz，这还只是理想情况。所以也能够得出一个结论：时域信号的采样率乘上一个固定系数即是变换之后频谱的宽度，可以用如下所示的一张图清晰说明：

采样率与频谱宽度的关系

经过上述的分析，我们有了如下的结论：更高的频谱分辨率需要有更长的采样时间，更宽的频谱分布需要提高对于原始信号的采样率，那我们在实际的使用过程中，当然是希望频谱更宽，分辨率更加精确，那么示波器的长存储就是必要的。

F103 如何进行 FFT

FFT 汇编库介绍

在本文的开头叙述了 ARM Cortex M4 具有 FPU 以及 DSP 指令，同时 ARM 官方也出了 DSP 方面的库来进行数字信号处理方面的工作，那么针对于 ARM Cortex M3 的 STM32F103 又是如何进行 FFT 的呢，显然，如果我们用 C 语言直接编写 FFT 算法，那样子的效率是极其低下的，因此，本文采用的方法是 ST 官方汇编 FFT 库的应用，由于官网现在找不到这个软件包，可以在公众号后台回复 FFT 获取软件包。
简单介绍一下，这个库是由汇编实现的，而且是基 4 算法，所以实现 FFT 在速度上较快。如果 X[N]是采样信号的话，使用 FFT 时必须满足如下两条：

N 得满足 4^n (n = 1,2,3…),也就是以 4 为基数。
采样信号必须是 32 位数据，高 16 位存实部，低 16 位存虚部(这个是针对大端模式)，小端模式是高位存虚部，低位存实部，一般常用的是小端模式。

汇编 FFT 的实现主要包括以下三个函数：

cr4_fft_64_stm32 : 实现 64 点 FFT
cr4_fft_256_stm32: 实现 256 点 FFT
cr4_fft_1024_stm32: 实现 1024 点 FFT

FFT 汇编库移植

将我们下载到的文件进行解压得到如下所示的文件：

解压得到的文件

进一步的我们需要将文件加入到我们的 keil 工程，加入工程之后的图如下所示：

汇编库添加

因为本文是针对于 256 点的 FFT ，因此只需要将cr4_fft_256_stm32 添加进来即可，加进来之后，再使用到 FFT 的文件里添加相关路径就可以。下面讲述具体的代码实例。

代码实例

FFT 计算幅值

首先我们定义采样的点数，以及 FFT 的输入数组，输出数组，以及各个谐波的幅值：

#define  NPT  256             /* 采样点数 */
uint32_t lBufInArray[NPT];    /* FFT 运算的输入数组 */
uint32_t lBufOutArray[NPT/2]; /* FFT 运算的输出数组 */
uint32_t lBufMagArray[NPT/2]; /* 各次谐波的幅值 */

在上述中，FFT 的输出数组和各次谐波的幅值的数组中只有采样点数的一半，是因为 FFT 计算出来的数据是对称的，因此通常而言只取一半的数据。
下面是波形采样并进行 FFT 的代码：

void adc_sample(void)
{
    for (i = 0; i < NPT; i++)
    {
        lBufInArray[i] = ADC_ConvertedValue[i];
    }
    cr4_fft_256_stm32(lBufOutArray,lBufInArray, NPT);
    GetPowerMag();
}

for循环里将 ADC 采集的数据存储到 FFT 运算的输入数组中去，在这里需要注意的是 STM32 是小端模式，因此采样数据是高位存虚部，低位存实部。紧接着就是调用汇编函数进行 FFT，FFT 运算之后，就进行幅值的计算，幅值的计算函数如下所示：

void GetPowerMag(void)
{
    signed short lX,lY;
    float X,Y,Mag;
    unsigned short i;
    for(i=0; i<NPT/2; i++)
    {
        lX  = (miniscope.freqency.lBufOutArray[i] << 16) >> 16;
        lY  = (miniscope.freqency.lBufOutArray[i] >> 16);

        //除以32768再乘65536是为了符合浮点数计算规律
        X = NPT * ((float)lX) / 32768;
        Y = NPT * ((float)lY) / 32768;
        Mag = sqrt(X * X + Y * Y)*1.0/ NPT;
        if(i == 0)    
            miniscope.freqency.lBufMagArray[i] = (unsigned long)(Mag * 32768);
        else
            miniscope.freqency.lBufMagArray[i] = (unsigned long)(Mag * 65536);
    }
}

上述代码中，lx 和 ly 的计算中，分别取的是FFT 的输出数组的高位和低位。进一步的，在计算 x 和 y 的时，除以 32768 是为了符合浮点数计算规律，至于为什么要进行浮点化，是因为浮点化就好像 10 进制里面的科学计数法。32768 = 2 的 15 次。除以 32768 也就是去除了浮点数后面的那个基数，只剩下前面的。比如 1991 改写成 1.991 * 10 的三次幂，除以 10 的三次方，只剩下 1.991，方便于下面的运算。而在后面又乘以 32768 和 65536 是因为要恢复到原先数据的大小，为什么下标为 0 的乘以 32768，而大于 0 的乘以 65535，是因为下标为 0 的代表的是直流分量，而剩余的是求出的乘以 2 才是实际模值。

FFT 计算频率

在本文的前面，笔者给出了这样一个公式用来计算 FFT 变换之后每个点对应的频率：
Fn = (n - 1) * Fs / N
N 是采样的点数，Fs 是采样频率，采样点数已经知道，还剩下采样频率未知，采样频率说白了也就是采样一个点的时间，也是 1 s 钟采样的点数，而这个该怎么确定呢。现有两种方法，第一种方法是在单片机进行 ADC 采集时，通过延时的方法每隔一段时间进行读取转换得到的数据，而这个延时的时间就是采样频率，这样听起来略显粗糙。另一种比较精确的方法，是通过 DMA + TIM 的方法，也就是通过 TIM 产生 PWM ，通过 PWM 触发 ADC 进行采集，这个时候，PWM 的频率也就是 ADC 的采样率，只需要控制 PWM 的频率就可以控制 ADC 的采样率，采集的数据通过 DMA 搬运至内存，当采样的点数达到规定的采样点数时，触发 DMA 中断，在中断里给出数据处理的信号，进一步进行 FFT，具体的原理及代码参考笔者的这篇文章：STM32 定时器触发 ADC 多通道采集，DMA搬运至内存。下面是一个简单的代码示例：

void wave_frequency_calculate(void)
{
    sample_frequency = 72000000.0 / (float)((sample_arr + 1) * (sample_psc + 1));
    wave_frequency = frequency_max_position * sample_frequency / NPT;
}

上述第一行代码是根据公式计算 pwm 的频率的公式，而 pwm 的频率也就是我们所需要的的采样频率。第二条代码中的 frequency_max_position 是除了直流分量幅值最大的点在数组中的位置，而这个点所对应的频率也就是我们采样波形的频率，至此，我们就计算出了采样波形的频率。