[TOC]

FFT函数的C语言实现和封装

一、序言

  本次封装脱离了ARM官方的DSP库（因为麻烦而且没有高点数的FFT）   我自己用C语言搓了一个FFT函数，为了实现更方便，空间复杂度更小的FFT。所以本驱动是适用于所有单片机型号的，唯一需要考虑的就是RAM的大小，RAM太小会导致数组开不大

二、定义说明

1.FPU

#define __FPU_PRESENT 1U    // 开启FPU

据说是开启FPU用的，可以提速，但根据实测好像并没有什么提升

2.NPT_max

#define NPT_max  16384

定义在FFTv2.h中，影响空间大小的主要因素。RAM不够了就把这个改小，但必须是2的整数次方

3.复数

typedef struct Complex
{
	float real; // 实部
	float imag; // 虚部
}Complex; // 复数

FFT要用到很多复数运算，为了便于组织，以结构体的形式存储复数

4.FFT运算结果

typedef struct FFTresult
{
	Complex result[NPT_max];    //FFT的结果，是复数形式
	float	Amp[NPT_max];       //幅值
	float	phase[NPT_max];     //相位
    float   dc;                 //直流分量
}FFTresult; //存放FFT运算结果的结构体

FFT的运算结果以复数的形式存储于result数组中 Amp存放result数组的模值 phase存放result数组的相位 dc存放的是直流分量（因为加窗之前要去掉直流，所以在这里记录下直流的大小）

5.函数

(1)初始化

void FFT_Init(void);

初始化函数，实际上只用来sin函数的打表，三角函数的计算很耗时间，查表可以节省很多时间（本驱动的运算时间与ARM官方的FFT运算时间相当）。在单片机上电之后调用一次就好

(2)FFT函数

/*
 *输入数据的类型为XXX的FFT
 *input: 输入数据
 *length: FFT点数，必须为2的整数次方
 *ch: FFTresult存放结果的结构体
 */
void FFT_u32(uint32_t* input, int length, FFTresult* ch);
void FFT_f32(float* input, int length, FFTresult* ch);
void FFT_int32(int* input, int length, FFTresult* ch);
void FFT_int16(int16_t* input, int length, FFTresult* ch);

FFT函数，为了兼容不同数据类型的输入（比如片上ADC是uint32_t类型，AD7606是int16_t类型），定义了4种函数（实际上这4个函数的内部代码是一样的，只不过我没有想到更好的兼容方法）

(3)窗函数生成

float Kaiser(int n, float beta, float* win);//Kaiser窗
float FlattopWin(int n, float* win);//Flattop窗
float Hanning(int n, float* win);//Hanning窗
float Hamming(int n, float* win);//Hamming窗
float rectangle(int n, float* win);//矩形窗

生成窗函数到win的地址，返回值是窗函数的恢复系数

(4)加窗FFT

/*
 *输入数据的类型为uint32_t的加窗FFT
 *input: 输入数据
 *length: FFT点数，必须为2的整数次方
 *win: 窗函数
 *ch: FFTresult存放结果的结构体
 */
void FFT_win_u32(uint32_t* input, int length, float* win, FFTresult* ch);
void FFT_win_f32(float* input, int length, float* win, FFTresult* ch);
void FFT_win_int32(int* input, int length, float* win, FFTresult* ch);
void FFT_win_int16(int16_t* input, int length, float* win, FFTresult* ch);

加窗FFT，win输入窗函数，无幅值恢复，需要额外调用幅值恢复函数

(5)幅值恢复

/*
 *恢复加窗后带来的幅值衰减
 *ch: 待恢复幅值的FFT结果
 *length: FFT点数
 *k: 恢复系数，生成窗函数时的返回值就是恢复系数
 */
void AmpRecovery(FFTresult* ch, int length, float k);

需要分析幅值的时候调用一下，如果只分析频率就不用管这个

(6)复数运算

Complex Complex_Multi(Complex c1, Complex c2);  // 复数乘法，返回值等于c1 * c2
Complex Complex_Add(Complex c1, Complex c2);    // 复数加法，返回值等于c1 + c2
Complex Complex_Minus(Complex c1, Complex c2);  // 复数减法，返回值等于c1 - c2
float Complex_Modular(Complex c);               // 复数求模
float Complex_arg(Complex c);                   // 复数求辐角

尽管可能用不到这些运算，为了保证灵活性，还是开放了这些函数的接口

三、调用步骤

1.不加窗的FFT

(1)引用头文件

#include "FFTv2.h"

(2)定义变量

uint32_t adc_buf[NPT_max];//这里的NPT_max也可以是比它小的常量
FFTresult ch1;//定义FFT结果变量

(3)初始化（打表提速）

FFT_Init();

(4)调用FFT

//根据你输入的数据类型选择不同的函数
FFT_u32(adc_buf, NPT, &ch1);

2.加窗FFT

(1)引用头文件

#include "FFTv2.h"

(2)定义变量

uint32_t adc_buf[NPT_max];//这里的NPT_max也可以是比它小的常量
float win[NPT_max];//定义窗函数，可扩展性体现在这里，你可以自己写窗函数，也可以用我提供的窗函数
FFTresult ch1;//定义FFT结果变量

(3)初始化（打表提速）

FFT_Init();

(4*)生成窗函数（不必要）

float k = FlattopWin(NPT, win);

(5)调用FFT

//根据你输入的数据类型选择不同的函数
FFT_win_u32(adc_buf, NPT, win, &ch1);

(6*)恢复幅值（不必要）

AmpRecovery(&ch1, NPT, k);

如果你不关注幅值的话，可以不调用这个

四、空间性能

主要占空间的是以下变量

typedef struct FFTresult
{
	Complex result[NPT_max];    //FFT的结果，是复数形式
	float	Amp[NPT_max];       //幅值
	float	phase[NPT_max];     //相位
    float   dc;                 //直流分量
}FFTresult; //存放FFT运算结果的结构体

uint32_t adc_buf[NPT_max];
float win[NPT_max];
FFTresult ch1;//定义FFT结果变量

float sin_table[NPT_max / 4];//sin打表

• result的大小：2 * NPT_max * 4 Byte

• Amp的大小：NPT_max * 4Byte

• phase的大小: NPT_max * 4Byte

• adc_buf的大小：NPT_max * 4 Byte

• win的大小：NPT_max * 4 Byte

• sin_table的大小：NPT_max Byte

不加窗的情况下占空间：21 * NPT_max Byte 加窗的情况下占空间：25 * NPT_max Byte