|
|
本帖最后由 nos002 于 2014-12-27 00:56 编辑
本人是个懒人，对于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[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;
- }
- }
- }
海涵原帖代码及完整代码出处：https://github.com/bjornfor/barcode_finder
|
本帖子中包含更多资源
您需要 登录 才可以下载或查看，没有帐号？注册
x
阿莫论坛20周年了！感谢大家的支持与爱护！！
月入3000的是反美的。收入3万是亲美的。收入30万是移民美国的。收入300万是取得绿卡后回国，教唆那些3000来反美的！
|
[复制链接]
发表于 2014-12-26 23:53:05
