中科因仑“3+1”工程特种兵精英论坛
标题:
实践“玩转FFT算法...任你移植”，正确AD采样及生成函数表
[打印本页]
作者:
张衍波
时间:
2015-10-7 12:53
标题:
实践“玩转FFT算法...任你移植”，正确AD采样及生成函数表
本人是个懒人，对于FFT算法，受水平限制自己是无论如何也理解不透，更写不出算法代码。好在有很多大神无私分享他们的成果，使得很多人站在大神的肩上（熟知的是巨人的肩上，好像不对，怎么能站到大神的肩上呢）忙碌自己的东西。
言归正传。最近手头有个东西需要显示交流电压和电流有效值，原本用的均方根法倒也简单，但硬件上没设过零检测无法自适应频率，于是搜了一通FFT测频，结果到了这里。虽然没找到测频的方法，但实践一下FFT测交流信号有效值也颇有兴趣。交流电压采样输入如下图：
输入波形不是很光滑，ADC为MCU的10位ADC，一周期采样64点。
搬运了海涵大神帖里的FFT代码，仔细阅读了相关注释，使用正弦函数生产软件生成了64点整数型正弦和余弦函数表。AD转换存入Fft_Real[]的数据为ADC的采样值减去直流电平值的绝对值。利用FFT计算的结果取基波Fft_Real[1]和Fft_Image[1]数据平方和后开根号，即基波的幅值，结果数据是随机的。生成64点全波正整形函数表，结果依然随机，百般折腾无果。没能力自己写出算法代码，只能充分信任搬运的了（╮(╯▽╰)╭）。
静下心来仔细看了几遍代码，发现FFT代码只取函数表的前半部的表值，尝试修改一些表达式以使能够取到全部表值，结果Fft_Real[]数据乱飞更凌乱了。分析整形正弦函数表前半部表值为正后半部表值为负，而代码只取前半部表值，算法思想应该是正、余弦波正负半波对称，那么要计算整个一周期的正、余弦函数，AD采样值就应该是一周波形的数据，而不是绝对值。据此修改了AD采样代码，正、余弦函数表采用32点半周整形数据。代码运行结果是基波幅值基本恒定。10位ADC量程700V，由于平移为单极性信号，使得ADC的分辨率降低一半，即只有700/512 V，代码运算结果取20次平均值后在正负1V跳动。只是FFT运算的幅值是实际的输入波形幅值的5倍，疑惑中，望高手解答！
此实践经验小结给同样的疑惑的电工，本人知识浅陋比较业余，高手见笑。
代码运行的硬件平台为PIC24FKL402，原样搬运的代码如下：
/**********************************************************************************************************************************************
* 这里讲讲傅立叶快速变换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<p24Fxxxx.h>
#include<math.h>
#include<stdlib.h>
#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
> Max)
Max = Vpeak
;
if (Vpeak
< Min)
Min = Vpeak
;
sum += Vpeak
;
}
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;
}
}
}
复制代码
海涵原帖代码及完整代码出处：
https://github.com/bjornfor/barcode_finder
转载
欢迎光临 中科因仑“3+1”工程特种兵精英论坛 (http://bbs.enlern.com/)
Powered by Discuz! X3.4