实践“玩转FFT算法...任你移植”，正确AD采样及生成函数表

胡小均 · 发表于 2015-5-16 10:53:32

本帖最后由胡小均于 2015-5-16 11:00 编辑

转载至阿莫论坛
做音乐流水灯的有福啦！！！

本人是个懒人，对于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，原样搬运的代码如下：

代码来了

[C] syntaxhighlighter_viewsource syntaxhighlighter_copycode

/**********************************************************************************************************************************************
* 这里讲讲傅立叶快速变换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[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;
        }
    }
}

胡小均 · 发表于 2015-5-16 11:16:39

自己顶！！！

李维强-15级 · 发表于 2015-5-16 19:33:22

有想法顶一个
可以设置为快速傅里叶变换，我想计算时间会短点

胡光鑫 · 发表于 2015-6-30 22:05:03

高端，顶一个

		自动登录	找回密码
密码			立即注册