适用于单片机的FFT快速傅里叶变换算法，51单片机都能用

• 源码
• • FFT.c
  • FFT.h
• 使用方法
• 效果
• 其他部分的代码
• • main.c

普中51-单核-A2
STC89C52
Keil uVision V5.29.0.0
PK51 Prof.Developers Kit Version:9.60.0.0

算法来自FFT算法的使用说明与C语言版实现源码 —— 原作者：吉帅虎

速度更快的版本见C语言实现的FFT与IFFT源代码，不依赖特定平台

移植十分简单，不依赖其他库，可自定义点数

#include 
#include "FFT.h"

struct compx Compx[FFT_N] = {0};		//FFT输入和输出：从Compx[0]开始存放，根据大小自己定义
double SIN_TAB[FFT_N / 4 + 1];			//定义正弦表的存放空间


struct compx EE(struct compx a, struct compx b)
{
	struct compx c;
	c.real = a.real*b.real - a.imag*b.imag;
	c.imag = a.real*b.imag + a.imag*b.real;
	return(c);
}


void create_sin_tab(double *sin_t)
{
	int i;
	for (i = 0; i <= FFT_N / 4; i++)
		sin_t[i] = sin(2 * PI*i / FFT_N);
}

double sin_tab(double pi)
{
	int n;
	double a = 0;
	n = (int)(pi*FFT_N / 2 / PI);

	if (n >= 0 && n <= FFT_N / 4)
		a = SIN_TAB[n];
	else if (n>FFT_N / 4 && n= FFT_N / 2 && n<3 * FFT_N / 4)
	{
		n -= FFT_N / 2;
		a = -SIN_TAB[n];
	}
	else if (n >= 3 * FFT_N / 4 && n<3 * FFT_N)
	{
		n = FFT_N - n;
		a = -SIN_TAB[n];
	}

	return a;
}

double cos_tab(double pi)
{
	double a, pi2;
	pi2 = pi + PI / 2;
	if (pi2>2 * PI)
		pi2 -= 2 * PI;
	a = sin_tab(pi2);
	return a;
}

void FFT(struct compx *xin)
{
	register int f, m, nv2, nm1, i, k, l, j = 0;
	struct compx u, w, t;

	nv2 = FFT_N / 2;					//变址运算，即把自然顺序变成倒位序，采用雷德算法
	nm1 = FFT_N - 1;
	for (i = 0; i < nm1; ++i)
	{
		if (i < j)						//如果i> (i != 0));
		xin[i].imag = i * sample_frequency / FFT_N;
	}
}


void Refresh_Data(struct compx *xin, int id, double wave_data)
{
	xin[id].real = wave_data;
	xin[id].imag = 0;
}

#ifndef FFT_H
#define FFT_H

#define FFT_N 16                                        //定义傅里叶变换的点数
#define PI 3.14159265358979323846264338327950288419717  //定义圆周率值

struct compx { double real, imag; };                    //定义一个复数结构

extern struct compx Compx[];							//FFT输入和输出：从Compx[0]开始存放，根据大小自己定义
extern double SIN_TAB[];								//正弦信号表

extern void Refresh_Data(struct compx *xin, int id, double wave_data);
extern void create_sin_tab(double *sin_t);
extern void FFT(struct compx *xin);
extern void Get_Result(struct compx *xin, double sample_frequency);

#endif

在FFT.h中修改 FFT_N 16，定义傅里叶变换的点数，由于FFT变换归一化后，除了0Hz的直流分量外，整个结果表是对称的，即若点数为16，则我们只看前8个点即可，所以这个傅里叶变换的点数可根据你的屏幕长方向上的像素数来决定，如128x64的屏幕可以考虑使用256点的FFT，我这里使用8x8的LED点阵屏来显示，故使用16点的FFT。

在运算前需调用一次 create_sin_tab(SIN_TAB) 建立正弦信号表, 之后将用查表法计算正弦值，加速计算

使用 Refresh_Data (Compx, i, wave_data)函数填入数据后

调用 FFT(Compx) 即可完成变换

使用 Get_Result (Compx, Sample_Frequency)函数求变换后结果的模值，存入复数的实部部分，频率存入复数的虚数部分，有效数据为前FFT_N/2个数

	#define Sample_Frequency 800			//采样频率 800Hz
	#define Frequency 100					//测试信号 100Hz

	for (i = 0; i < FFT_N; ++i)				//使用Refresh_Data(Compx, i, wave_data)函数填入数据，这里建立一个100Hz的方波测试信号
	{ 
		if (sin(2 * PI * Frequency * i / Sample_Frequency) > 0)
			Refresh_Data(Compx, i, 1);
		else if (sin(2 * PI * Frequency * i / Sample_Frequency) < 0)
			Refresh_Data(Compx, i, -1);
		else
			Refresh_Data(Compx, i, 0);
	}										

	create_sin_tab(SIN_TAB);				//建立正弦信号表, 之后将用查表法计算正弦值，加速计算
	FFT(Compx);								//快速傅里叶变换
	Get_Result(Compx, Sample_Frequency);	//求变换后结果的模值，存入复数的实部部分，频率存入复数的虚数部分，有效数据为前FFT_N/2个数

计算频率的方法：
采样频率为800Hz，共16个点，故一格 = 800Hz / 16 = 50 Hz
如图，在第三格和第七格各有一个峰，则对应的频率是2 x 50Hz = 100Hz 和 6 x 50Hz = 300Hz，该结果已由Get_Result函数计算并存在原数组的虚部中。

物理意义的解读见如何 FFT(快速傅里叶变换) 求幅度、频率（超详细 含推导过程）—— Xav Pan

性能
晶振频率 11.0592MHz 6T模式（22.1184MHz 12T）
仿真下花费约 0.13222819 - 0.06633789 = 0.0658903 (s)
即完成一次16点的FFT变换，11.0592MHz 6T的51单片机花费 65.8903 ms


同样的程序，使用Vs 2015 跑65536个点的计算结果：

同样的数据使用python计算800个点的结果：
来源：使用python（scipy和numpy）实现快速傅里叶变换（FFT）最详细教程 —— LoveMIss-Y

点阵屏部分代码见【51单片机快速入门指南】2.4：IO扩展(串转并) 与 LED点阵屏

#include 
#include "intrins.h"
#include "stdint.h"
#include "FFT.h"
#include 
#include "HC74595.h"

void main(void)
{
	uint8_t i, j; 
	double Largest = 0;

	#define Sample_Frequency 800			//采样频率 800Hz
	#define Frequency 100					//测试信号 100Hz

	for (i = 0; i < FFT_N; ++i)				//使用Refresh_Data(Compx, i, wave_data)函数填入数据，这里建立一个100Hz的方波信号
	{ 
		if (sin(2 * PI * Frequency * i / Sample_Frequency) > 0)
			Refresh_Data(Compx, i, 1);
		else if (sin(2 * PI * Frequency * i / Sample_Frequency) < 0)
			Refresh_Data(Compx, i, -1);
		else
			Refresh_Data(Compx, i, 0);
	}										

	create_sin_tab(SIN_TAB);				//建立正弦信号表, 之后将用查表法计算正弦值，加速计算
	FFT(Compx);								//快速傅里叶变换
	Get_Result(Compx, Sample_Frequency);	//求变换后结果的模值，存入复数的实部部分，频率存入复数的虚数部分，有效数据为前FFT_N/2个数

	for (i = 0; i < FFT_N / 2; ++i)			//将计算结果处理成图像的部分
	{         
		if(Compx[i].real > Largest)
			Largest = Compx[i].real;
	}
	for(i = 0; i < 8; ++i)
	{
		for(j = 0; j < (uint8_t)((Compx[i].real / Largest) * 8 + 0.5); ++j)
			Mat[8 - j - 1][i] = 1;
	}

	while(1)
	{
		imshow(Mat);						//显示图像
	}
}