/**********************************************************************************************************************************************

* 这里讲讲傅立叶快速变换FFT算法的C例程，不恰当的地方还请探讨。

* 此例程可以移植到任何的单片机，但前提是单片机内存要够大(如果做32点的应该256字节内存就足够了)，至于速度快或慢，只要你能承受那些慢速的就无所谓。

* 这里以128点采样作为例子。要注意的是，采样点数N和采样率Fs还有信号频率F的关系。

* FFT算法可以把对波形采样回来的数据进行转换，比如有一个模拟信号过来，假设对其以某个采样率采样了128个点的数据，经过FFT换算以后，就可以获得这段

* 波形里面含有的不同频率的特征，也就是变成频谱。换算的结果是以复数的形式表现的，比如某个频率的信息以a＋bi的形式表示，这里a是实部，b是虚部。

* 而例程中会以数组s16 Fft_Real[128]来表示128个频率的实部，以s16 Fft_Image[128]来表示128个频率的虚部。

* 采样点数N和采样率Fs还有信号频率F的关系是这样的，在一个模拟信号中，采样率至少应该是这个信号中需要分析的最高频率的2倍，而这个信号中需要分析的

* 最低频率Fl＝Fs/N，也就是采样率除以采样点数，同时这个频率也是fft换算过程的频率分辩率。可能有点晕了，这里举例说明。

* 假设我们有一个信号里面含有100Hz的信号和含有1000Hz的信号叠加在一起，用示波器是很难分清的，通过对这信号进行采样和FFT转换就可以得到这组信号的

* 频谱特征，因为需要分析的最低信号为100Hz，我们假设现在对信号要进行128点的采样，那么我们至少要能分辨100Hz的信号，采样率可以是100Hzx128点＝12800Hz，

* 也就是我们需要以每秒采集12800点的速度来采集1/100秒的时间(因为只采集128点，所以时间上是1/100秒)，当然也可以使用6400Hz的采样率，这样分辩率是50Hz，

* 这里就以12800Hz为例。

* 采集回来的数据要怎么处理呢？首先对于FFT算法，需要有数组数据输入，经过FFT换算以后会以数组数据输出。我们可以直接定义两个数组s16 Fft_Real[128]和

* s16 Fft_Image[128]，我们把采集来的128个数据按一个特殊的顺序输入到s16 Fft_Real[128]这个数组里面(当作实部输入)，而s16 Fft_Image[128]则全部输入0，

* 然后运行fft算法转换。转换过后，s16 Fft_Real[128]里面的128个数会变成128个频率的实部，而s16 Fft_Image[128]里面的128个数会变成128个频率的虚部。

* 哪128个频率呢，实际上只能获得64个频率。在这两个输出的数组里面，数组s16 Fft_Real[128]的第一个数(第一个数是s16 Fft_Real[0])和s16 Fft_Image[128]

* 的第一个数(第一个数是s16 Fft_Image[0])分别是0Hz这个频率的复数形式结果里面的实部和虚部，而第二个数是100Hz的复数形式的结果，第三个数是200Hz的结果，

* 依此类推，第64个数是6300Hz的结果。假如我们想知道信号中是否含有100Hz这个频率存在，那么我们需要查看s16 Fft_Real[1]和s16 Fft_Image[1]这两个数，

* 把实部和虚部分别进行平方求和以后再开平方，就是100Hz信号的模，把模值除以采样点数N的1/2就得到100Hz的信号的峰值。同理，要知道信号里面是否含有1000Hz，

* 我们需要取出s16 Fft_Real[10]和s16 Fft_Image[10]来计算。这128组数据里面，只有前面64组是有用的，后面的都是镜像数据，我们不理会。因为采样率是12800Hz，

* 我们能查看的最高频率是6400Hz，也就是需要拿s16 Fft_Real[64]和s16 Fft_Image[64]出来算，但这两个值很可能不再准确，因为是极限的问题。

*

* 以下我们来看看相关的例程，这个例程可以用于移植到任何的单片机中，只要单片机的速度不太慢(太慢你愿意慢慢等也没关系)，内存足够就可以运行。这个例程可以

* 修改为256点和512点或者1024点甚至更高点数的fft，下面会讲到需要修改一些什么参数。增加采样点数可以在采样率不变的情况下增加分辩率，如果点数增加一倍，

* 采样率也提高一倍，那么分辩率不变，但能识别的最高频率将提高一倍，这是采样点数增加的好处，但同时会减慢运算速度。此例程使用stm32f103在72M主频下测试256点fft,

* 运算时间为10mS左右，官方的函数库不知道会不会经过优化而更快，但考虑到许多人不喜欢使用函数库，也考虑这个例程的可移植性，这里直接使用fft函数而不讲函数库的内容。

*********************************************************************************************************************************************/

#include

#include

#include

#include"GenericTypeDefs.h"

#define NUM_FFT 64  //这里要算多少点的fft就赋值多少，值只能是2的N次方

#define N       6   //这里因为128是2的7次方，如果是计算256点，则是2的8次方，N就是8，如果是512点则N＝9，如此类推

//以下为FFT傅立叶转换算法用到的相关定义

//UINT8 ADC_Count = 0;

/********注意，这里u8就是8位的unsigned char类型数据，在51里面就是unsigned char。以下的数组s16 Fft_Real[128]就是16位signed类型的数组，相当于signed int，

* 而Fft_Real只是一个数组的名字，随便起也可以的，Fft_Image也一样**********/

//这是用到的一些寄存器，注意这里如果要算256点以上的fft的话，需改成16位的u16。

//中间临时变量，名称也是自己定义的，但要与fft函数里面的对应

//UINT16 TEMP1; //用于求功率的，可不需要

INT16 Fft_Real[NUM_FFT]; //fft实部，128数组

INT16 Fft_Image[NUM_FFT]; //fft虚部，128数组

UINT8 FFTSwitch = 0;

UINT8 CalculateSW = 0;

float Vrms;

void Fft_Imagclear(void);

/*********************************************************************************************************************************

* 以下为放大128倍后的sin正弦函数数组表格，这个表格可以用“正弦波表生成”这个软件来生成，要注意需要做多少点的fft，就需要生成多少点的。

* 数据类型要选择整形，不要选择正整型，这里为了能在普通单片机快速运行，不使用浮点数，使用查表法以达到最快的运算。如果是256点以上的，表格

* 数据类型s8要改成s16的，以下的两个数组表格也同样如此，注意这里相当于指针用法，如果内存足够的，最好把表格写入内存，那样运行速度快，如果

* 内存资源紧的，就把表格写入程序区，对于51就是一个signed char code table[N]和signed char table[N]的区别，带了code就告诉编译器我这个表

* 格要放在ROM程序存储区，不带code就直接放入内存，其他单片机不同写法自己研究研究

*

* *******************************************************************************************************************************/

const INT8 SIN_TAB[32] = {0x00, 0x0c, 0x18, 0x24, 0x30, 0x3b, 0x46, 0x50,

0x59, 0x62, 0x69, 0x70, 0x75, 0x79, 0x7c, 0x7e,

0x7f, 0x7e, 0x7c, 0x79, 0x75, 0x70, 0x69, 0x62,

0x59, 0x50, 0x46, 0x3b, 0x30, 0x24, 0x18, 0x0c};

//以下是放大128倍后的cos余弦函数数组表格，这里注意事项与上面相同，只不过选择余弦来生成

const INT8 COS_TAB[32] = {0x7f, 0x7e, 0x7c, 0x79, 0x75, 0x70, 0x69, 0x62,

0x59, 0x50, 0x46, 0x3b, 0x30, 0x24, 0x18, 0x0c,

0x00, -0x0c, -0x18, -0x24, -0x30, -0x3b, -0x46, -0x50,

-0x59, -0x62, -0x69, -0x70, -0x75, -0x79, -0x7c, -0x7e};

/************************************************************************************************************

* 以下是采样存储序列表，这个表格可以在FFT_Code_Tables.h这个文件里面找到，更多点的FFT也在里面找，就是里面的unsigned

* int code BRTable[N]的那些，是用来控制采样点数据输入到实部数组s16 Fft_Real[128]里面的

************************************************************************************************************/

const UINT8 LIST_TAB[64] = {

0, 32, 16, 48, 8, 40, 24, 56,

4, 36, 20, 52, 12, 44, 28, 60,

2, 34, 18, 50, 10, 42, 26, 58,

6, 38, 22, 54, 14, 46, 30, 62,

1, 33, 17, 49, 9, 41, 25, 57,

5, 37, 21, 53, 13, 45, 29, 61,

3, 35, 19, 51, 11, 43, 27, 59,

7, 39, 23, 55, 15, 47, 31, 63

};

/********************************************************************************************************************************

* 以下为fft算法主体部分

* ******************************************************************************************************************************/

void FFT(void)

{

UINT8 i, j, k, b, p;

INT16 Temp_Real, Temp_Imag, temp;

UINT16 Max, Min;

UINT32 sum;

static UINT16 Count = 0;

static UINT16 Vpeak[30];

if (FFTSwitch == 1)

{

Fft_Imagclear();

for (i = 1; i <= N; i++) /* for(1) */

{

b = 1;

b <<= (i - 1); //蝶式运算，用于计算 隔多少行计算。例如第一级 1和2行计算，，，第二级

for (j = 0; j <= b - 1; j++) /* for (2) */

{

p = 1;

p <<= (N - i);

p = p*j;

for (k = j; k < NUM_FFT; k = k + 2 * b) /* for (3) 基二fft */

{

Temp_Real = Fft_Real[k];

Temp_Imag = Fft_Image[k];

temp = Fft_Real[k + b];

Fft_Real[k] = Fft_Real[k] + ((Fft_Real[k + b] * COS_TAB[p]) >> 7) + ((Fft_Image[k + b] * SIN_TAB[p]) >> 7);

Fft_Image[k] = Fft_Image[k] - ((Fft_Real[k + b] * SIN_TAB[p]) >> 7) + ((Fft_Image[k + b] * COS_TAB[p]) >> 7);

Fft_Real[k + b] = Temp_Real - ((Fft_Real[k + b] * COS_TAB[p]) >> 7) - ((Fft_Image[k + b] * SIN_TAB[p]) >> 7);

Fft_Image[k + b] = Temp_Imag + ((temp * SIN_TAB[p]) >> 7) - ((Fft_Image[k + b] * COS_TAB[p]) >> 7);

//移位，防止溢出。结果已经是本值的1/64

Fft_Real[k] >>= 1;

Fft_Image[k] >>= 1;

Fft_Real[k + b] >>= 1;

Fft_Image[k + b] >>= 1;

}

}

}

Vpeak[Count++] = (UINT16)(100 * sqrt(Fft_Real[1] * Fft_Real[1] + Fft_Image[1] * Fft_Image[1]));

if (Count == 30)

{

sum = 0;

Max = Vpeak[0];

Min = Vpeak[0];

for (i = 0; i < Count; i++)

{

if (Vpeak[i] > Max)

Max = Vpeak[i];

if (Vpeak[i] < Min)

Min = Vpeak[i];

sum += Vpeak[i];

}

Vrms = (sum - Max - Min)/ (Count - 2);

Vrms = Vrms / 32 /1.414;

Count = 0;

}

FFTSwitch = 0;

}

//         Fft_Real[0]=Fft_Image[0]=0;  //去掉直流分量，也可以不去掉

//   Fft_Real[63]=Fft_Image[63]=0;

//以上已经把128点的实部和虚部求完，下一次运算前需要把所有虚部重新清零。

//要求某个频率点的模，则模值=根号(实部平方+虚部平方)，即sqrt((Fft_Real[n]*Fft_Real[n])+(Fft_Image[n]*Fft_Image[n]))，

//第n个频率点的值是数组上的Fft_Real[n]和Fft_Image[n]，而Fft_Real[0]是直流分量。Fft_Real[1]是最低频率点，也是最小频率分辩率值，等于采样率/采样点数N。

//波形峰值大小=模值/(N/2)， N为采样点数。

//注意，由于上面把求得的值已经移位除法，相当于缩小了64倍(移位7位好像是128倍吧？？后面为什么还要移动一位？这里待高手指点，本人也不是很清楚，这里只做移植总结)。

//得到的那些实部虚部的结果爱怎么处理怎么处理，可以做音频频谱强度显示啦什么的。

}

/************************************************************************************************************************************

*

************************************************************************************************************************************/

void Fft_Imagclear(void) //fft虚部清零函数，在运行FFT函数之前需要先运行这个

{

volatile UINT8 cnt;

for (cnt = 0; cnt < NUM_FFT; cnt++)

{

Fft_Image[cnt] = 0;

}

}

// 以下是和MCU硬件相关的部分

/********************************************************************************************

* 模 块：ADC初始化

********************************************************************************************/

void ADC_Init(void)

{

ANSA |= 0x0002; // RA1为模拟输入

TRISA |= 0x0002;

AD1CON1 = 0x0000;

AD1CON2 = 0x0000;

AD1CON3 = 0x0000;

AD1CHS = 0x0000;

AD1CSSL = 0x0000;

AD1CON1bits.SSRC = 0b010; // 由Timer1 比较结束采样并启动转换

AD1CON1bits.ASAM = 1; // 1 = 最后一次转换结束后立即开始采样； SAMP 位自动置1

AD1CON1bits.ADSIDL = 1;

AD1CON2bits.VCFG = 0b000; // VREF+接AVDD,VREF-接AVSS

AD1CON2bits.OFFCAL = 0; // 获取实际值

AD1CON2bits.CSCNA = 0; // 不扫描输入

AD1CON2bits.SMPI = 0b0000; // 每完成1个采样产生中断

AD1CON2bits.ALTS = 0; // 总是使用MUX A输入多路开关设置

AD1CON3bits.ADCS = 0b00001; // 转换时钟 2TCY

AD1CON3bits.SAMC = 0b00100; // 自动采样时间 8*TAD

AD1CHSbits.CH0SA = 1; // 选择AN1

AD1CHSbits.CH0NA = 0;

IPC3bits.AD1IP = 0b011;

IEC0bits.AD1IE = 1;

IFS0bits.AD1IF = 0;

AD1CON1bits.ADON = 1; //打开AD模块

}

/**************************************************************

* 模 块：TIMER1初始化

***************************************************************/

void T1_Init(void)

{

TMR1 = 0;

PR1 = (UINT16)((16129.0/SAMPLE_POINTS)/Tcy); //5000 = (20ms*1000/64)/0.0625

T1CONbits.TCS = 0; //内部时钟，Fosc/2

T1CONbits.TCKPS = 0b00; //预分频 1:1

T1CONbits.TGATE = 0; //0 = 禁止门控时间累加

IPC0bits.T1IP = 0b011;

IFS0bits.T1IF = 0;

//IEC0bits.T1IE = 1;

T1CONbits.TON = 1;

}

/*************************************************************************************

* 模 块：TIMER1

* 用 途：用于结束ADC采样并启动转换，由ADC1CON1设置

*************************************************************************************/

void __attribute__((interrupt, shadow, no_auto_psv)) _T1Interrupt()

{

IFS0bits.T1IF = 0; // manually cleared MI2C1 Interrupt flag

}

/*************************************************************************************

* 模 块：ADC中断

* 用 途：每中断一次读取一个采样点的瞬时值

*************************************************************************************/

extern UINT8 FFTSwitch;

void __attribute__((interrupt, shadow, no_auto_psv)) _ADC1Interrupt(void)

{

static UINT8 i = 0,j = 0;

IFS0bits.AD1IF = 0;

SampleBuf[i++] = ADC1BUF0;

if (SAMPLE_POINTS == i)

{

i = 0;

SampleCycleMsg = 1;

}

if(FFTSwitch == 0)

{

Fft_Real[LIST_TAB[j++]] = ADC1BUF0 - MEDIAN_VREF; // ADC采样值 - 参考电压半值，即得到交流采样信号的正负半波瞬时值

if (SAMPLE_POINTS == j)

{

j = 0;

FFTSwitch = 1;

}

}

}