2023电赛H题详解1 - FFT算法

UID196129582

已于 2024-04-01 04:53:37 修改

阅读量2.9k

点赞数 13

分类专栏：天天向上的电赛路文章标签：算法 stm32 c语言单片机傅里叶分析

于 2023-08-15 07:38:03 首次发布

本文链接：https://blog.csdn.net/UID196129582/article/details/132289005

版权

天天向上的电赛路专栏收录该内容

2 篇文章 0 订阅

订阅专栏

FFT，即快速傅里叶变换，是一种测量信号频率的方法。在2023年H题、2021年A题中，均用到了fft算法。FFT算法往往配合高速AD采样、DMA配合使用，由于stm32cubemx里面有DSP库，所以我们可以使用官方的fft函数，这样做简单高效，准确度高。

下面我来详细介绍一下FFT算法，以及如何使用ADC+FFT测量信号频率。我争取用这一篇文章，帮你们彻底学会stm32的FFT实现。

我的2023H题代码在这里，建议看懂这篇文章，再根据我的代码修改。直接用大概率会出问题，毕竟我们用的开发板可能不一样，而且我用了三个AD采样引脚（实际只用了两个，第三个本来打算做反馈调节的，时间不够就没搞了），由接了3个AD9833模块作为输出，所以直接复制粘贴肯定搞不定。

链接: https://pan.baidu.com/s/1dMlcS2DNF-SYfLU5yADhtQ?pwd=ee7n 提取码: ee7n 复制这段内容后打开百度网盘手机App，操作更方便哦

FFT详解

什么是FFT

FFT是快速傅里叶变换的简称，这里我不介绍数学上的原理了，想要了解数学原理的话可以看这个：十分简明易懂的FFT（快速傅里叶变换）_路人黑的纸巾的博客-CSDN博客

咱们搞单片机的可以不需要知道它是怎么算的，只需要知道FFT的作用就行——计算信号的频率。对傅里叶变换有了解的都知道任何函数都可以写成几个正弦函数的和，只要知道各个sin函数的模值、频率、相位，就可以将这些个正弦函数写出来，把他们加起来就是这个信号的波形。

FFT可以得到不同频率的sin对应的实部和虚部，根据实部和虚部就可以得到模值与角度，从而写出信号的傅里叶变换式。咱们学单片机的只需要知道怎么处理FFT得到的数据，从而计算出信号频率就OK了。

单片机FFT原理及使用

单片机中，FFT通常会与AD采样一起使用。它的原理很简单，如下图：

首先要确定采样频率和采样个数。根据奈奎斯特采样定理，采样频率需要大于信号频率的两倍。由于采样频率越低，采集相同的点所需的时间越长，采集到的周期数会越多，结论会越精确，所以理论上采样频率越低越好，且不能低于输入信号的两倍。

采样个数的话，一般我是取1024个点。理论上采样的点越多越精确，越少花的时间越短。不过采样点数太少会导致LSB特别大，精度降得很厉害。采样点数太多的话会导致采样时间过长，且占用的内存会很大，计算量也会比较大。采样点数是2的k次幂，官方函数貌似最大是4096个点（我没试过）。

确定了采样频率f_s和采样个数n后，就可以使用ADC采样了。由于采样的点较多，建议使用DMA传输ADC数据。将这n个数据放入到数组a[n]里，这个数组a[n]就是输入的数据。然后我们就可以开始进行FFT计算了。

需要注意的是，官方FFT函数需要输入的数组是A[2n]，这需要你对a[n]进行预处理。将ADC数据按顺序放到A[2n]的偶数位作为实部，即 A[2k]=a[k] ；将0放到A[2n]的奇数位作为虚部，即 A[2k+1]=0 。

使用官方的FFT函数arm_cfft_f32，可以得到数组b[2n]。其中的偶数位对应的是实部，奇数位对应的是虚部，n对应的是频率。也就是说b[2n]其实就是讲信号进行了傅里叶变换，得到了不同频率的正弦函数，只不过他们是用实部和虚部表示的。

大部分情况下，我们不需要知道相位，只需要知道模值就行。于是我们可以使用官方的FFT函数arm_cmplx_mag_f32，得到对应的模值c[n]。说白了就是将n对应的实部和虚部开平方相加再开根号。

总而言之，单片机的FFT算法干了这样一件事，将ADC采到的n个数进行了傅里叶变换，得到了

$f(t)=a_{0} + a_{1}\sin \omega _{1}t + a_{2}\sin \omega _{2}t + ... + a_{n}\sin \omega _{n}t$

其中a0、a1、a2就是对应频率的模值，数组c[n]中的n与 $\omega _{n}$ 一一对应。由于AD采样必然存在误差，所以每一个 $a_{x}$ 都不为0，但是有大有小。如果输入的信号是10kHz的正弦波，那么10kHz对应的x值所对应的 $a_{x}$ 就会非常大，远大于别的模值。 $a_{x}$ 的大小与信号的幅度有关，幅度越大， $a_{x}$ 越大。 $a_{0}$ 则对应的是信号直流分量的模值。

如果是1Vpp、50kHz的信号A 与 1Vpp、100kHz的信号B 相加得到的信号C，进行ADC+FFT后，50kHz与100kHz对应的 $a_{x_{1}}$ 与 $a_{x_{2}}$ 就会很大，其他的 $a_{x}$ 就会很小。

FFT函数

注意：使用FFT函数需要导入DSP库，CLion如何操作可以看我后面的实操。Keil的操作方法可以在CSDN上搜。

arm_cfft_f32

该函数可以这样使用：

arm_cfft_f32(&arm_cfft_sR_f32_len1024, fft_inputbuf, 0, 1);

第一项表示FFT类型，取1024个点的FFT算法就用&arm_cfft_sR_f32_len1024。需要注意的是，CLion里不能直接用&arm_cfft_sR_f32_len1024，需要在main函数前面加上：

extern const arm_cfft_instance_f32 arm_cfft_sR_f32_len1024;

第二项fft_inputbuf就是输入的数组A[2n]。注意，A[2n]是a[n]经过处理得到的（奇数位置0，偶数位依次填入），且须定义为float型。

第三项和第四项就填0,1即可。

该函数直接对A[2n]进行处理，即运行该函数后，fft_inputbuf即为输出数组b[2n]，偶数位表示实部，奇数位表示虚部。

这个函数有点复杂，我找到了一个比较详细的介绍，由于当时只截了一个图，所以没法注明出处了。如果原作者看到了可以联系我。

arm_cmplx_mag_f32

该函数可以这样使用：

arm_cmplx_mag_f32(fft_inputbuf, fft_outputbuf, FFT_LENGTH);

第一项表示输入数组b[2n]，第二项表示输出数组c[n]，第三项表示采样点数（1024）。其实就是把实部和虚部的平方和开根号。

以上就是理论部分，下面是使用 stm32cubeMX + CLion 对 ADC+FFT 的实操。我个人比较推荐有C语言基础的同学使用CLion编程，因为Keil的自动补全功能不好用，且调试起来比较复杂。至于如何搭建CLion，可以参考：从零开始教你使用Clion优雅开发STM32（一）软件安装与环境配置_clion开发单片机_王拉图的博客-CSDN博客

当然，如果你在使用CLion的过程中遇到了问题，欢迎在评论区提问。

FFT实操

其实已经过去好几个月了，我几乎都忘记这个没更完了。不过居然真的有同志催我更新，那我先更着，争取这几天更完，不然真忘干净了。

ADC配置

说起来ADC的配置啥的，时间真的过去太久了，我不敢误人子弟，实操步骤就把我打电赛的时候看的教程给你们看吧。

首先，要弄明白ADC是啥？很简单，A是模拟信号，C是数字信号，ADC就是单片机里实现将数字信号转化为模拟信号的功能。

模拟信号就是单片机常用的0~3.3V的电压，单片机发出的方波、正弦波之类的都是模拟信号，是不是很好懂？数字信号是具体的数值，在stm32里就是0~4095这些数，这些数字与0~3.3V的电压是对应的。

具体的计算方式为：

ADC值 = （实际电压值 / 3.3）* 4095 。当然，由于单片机识别电压会存在一定的误差，所以不一定100%准确。

弄清楚了这个，配置啥的都简单了。看看我推荐的这两个教程就OK了。

教程1：STM32CUBEMX配置教程（十二）STM32的定时器触发的固定频率ADC采样（使用DMA）_定时器溢出踩adc-CSDN博客

这个教程是从头到尾教你ADC是什么、怎么弄的，很详细，推荐小白学一学，基本上花1h就能彻底搞明白。

当然，你要是真没时间，咱也理解，所以贴心的为你找了ADC+FFT+DAC的全套流程。

STM32 ADC+定时器+DMA+FFT_stm32adc采样 fft-CSDN博客

我当时就是跟着这个学的，需要注意的是，一定要在ADC回调函数中加入回调标志！

另外，他这个教程是使用keil的，但这年头谁还用keil啊，如果你跟我一样用的是clion，那么一定要记得导入DSP库，不然无法使用fft函数。

DSP库导入

导入DSP库我就直接放别人的攻略了，自己跟着弄一遍就行了。

在Clion中配置DSP库_clion dsp-CSDN博客

不用担心这个教程有问题，我就是看这个教程弄的。出了问题一定是你自己的原因，实在解决不了可以在评论区里问。

FFT代码

main函数的代码：

FFT_deal0(ad_value, fft_in_adc_data_1, fft_adc_n);
arm_cfft_f32(&arm_cfft_sR_f32_len1024,
             fft_in_adc_data_1, 0, 1);
arm_cmplx_mag_f32(fft_in_adc_data_1, fft_out_adc_data_1, fft_adc_n);

各个变量的含义：（main函数的定义和前期处理）

#define fft_adc_n 1024    // 这里采用的1024点FFT，也可以换成256之类的

uint16_t ad_value[fft_adc_n * 2][3] = {0};
// ad_value是你前面ADC采样得到的值
// ADC采样频率需要大于奈奎斯特采样定理的f_min，理论上讲f_s越贴近f_min，FFT精确度更高


float fft_in_adc_data_1[fft_adc_n * 2] = {0}; // FFT输入，实部是ADC值，虚部都是0
// 这个值是FFT_deal0的返回值，不用你干什么


// &arm_cfft_sR_f32_len1024 是DSP库里关于fft的依据，你有这个就能实现1024点的fft


float fft_out_adc_data_1[fft_adc_n] = {0};  // FFT输出
// 这个是处理完后返回的fft值，可以直接打印看看结果

FFT_deal0函数：

extern uint16_t ad_value[fft_adc_n * 2][3];
extern uint16_t id[fft_adc_n];

void FFT_deal0(uint16_t ad_value[][3], float *fft_in_1, int data_length) {
    int i;
    for ( i = 0; i < data_length; i++ ) {
        fft_in_1[2 * i] = (float) ad_value[i][0];
        fft_in_1[2 * i + 1] = 0;
    }
}

冒泡排序（后面计算FFT极大值时要用，FFT极大值表示的就是对应波的频率）

void BubbleSort(float k[], int n, uint16_t id[]) {
    int i, j, flag = 1;
    float temp;
    uint16_t tmp;
    for ( i = n - 1; i > 0 && flag == 1; i-- ) {
        flag = 0;
        for ( j = 1; j < i; j++ )
            if ( k[j] >= k[j + 1] ) {
                temp = k[j];
                k[j] = k[j + 1];
                k[j + 1] = temp;
                tmp = id[j];
                id[j] = id[j + 1];
                id[j + 1] = tmp;
                flag = 1;
            }
    }
}

后续处理参考（主要是为了验证FFT的准确性，放main函数fft代码后面进行检验）

uint16_t id[fft_adc_n];

for ( i = 0; i < fft_adc_n; i++ ) {
    printf("%d %.2f %.2f\n", i, fft_in_adc_data_1[2 * i], fft_in_adc_data_1[2 * i + 1]);
}
printf("xxxxxxxxxxxxxxxx\n");
printf("xxxxxxxxxxxxxxxxx\r\n");
for ( i = 0; i < fft_adc_n; i++ ) {
    id[i] = i;
}
BubbleSort(fft_out_adc_data_1, fft_adc_n / 2, id);

后续处理（这是我个人的处理，仅供参考）

float f_a，f_b;      //  A和B的频率
float A_a，A_b;      //  A和B的幅度


    int f_pre = 0;
    for ( i = 0; i < fft_adc_n / 2; i++ ) {
        if ( id[i] < fft_adc_n / 2 ) {
            printf("%d  %.2f\r\n", id[i], fft_out_adc_data_1[i]);
        }
    }
    if ( id[i - 1] < id[i - 2] ) {     //  找到�??????????高频�??????????
        f_pre = id[i - 2];
        f_a = (float) id[i - 1] / fft_adc_n * f_s;
        A_a = fft_out_adc_data_1[i - 1];
        f_b = (float) f_pre / fft_adc_n * f_s;
        A_b = fft_out_adc_data_1[i - 2];
    } else {
        f_pre = id[i - 1];
        f_a = (float) id[i - 2] / fft_adc_n * f_s;
        A_a = fft_out_adc_data_1[i - 2];
        f_b = (float) f_pre / fft_adc_n * f_s;
        A_b = fft_out_adc_data_1[i - 1];
    }