打印
[STM32F1]

stm32做FFT变换

[复制链接]
10674|11
手机看帖
扫描二维码
随时随地手机跟帖
跳转到指定楼层
楼主
892953881|  楼主 | 2015-12-29 21:02 | 只看该作者 |只看大图 回帖奖励 |倒序浏览 |阅读模式
用AD采样后对电压进行FFT变换。在网上找了很多算法无非都是这样的差不多的。假如用STM3210ZET6 12位AD采样  1024个点,采样频率1024HZ.
这个函数的输入是个结构体表示的是一个复数,而我做AD采样后的是一个12位的值,请问下各位如何将AD采样后的值转换成这个函数输入所需要的。或者哪位有其他算法的代码求分享下。谢谢了!!!


compx fftres[FFT_N]; //FFT数据段                                   
//m^n函数
u32 mypow(u8 m,u8 n)
{
        u32 result=1;         
        while(n--)result*=m;   
        return result;
}   
//快速傅里叶变换
//32/64/128/256/512/1024点的FFT
//STM32 计算1024点费时35.7ms左右@72M
//如果超频到120M,则时间只需要22ms左右了
//N:傅里叶变换的点数
//xin:输入数组大小为N+1
void FFT(compx *xin,u16 N)
{
        u16 f,m,LH,nm,i,k,j,L;
        u16 p;
        u16 le,B,ip;                      
        compx ws,t;
        LH=N/2;
        f=N;
        for(m=1;(f=f/2)!=1;m++);//求得M的值
        nm=N-2;   
        j=N/2;                  
        for(i=1;i<=nm;i++)//码位倒置
        {
                if(i<j){t=xin[j];xin[j]=xin;xin=t;}//码位家换
                k=LH;
                while(j>=k){j=j-k;k=k/2;}
                j=j+k;
        }   
        for(L=1;L<=m;L++) //fft  傅里叶变换
        {  
                le=mypow(2,L);//用自己定义的乘方函数,效率比库函数高很多.这里如果采用移位计算,效率更高.          
                B=le/2;                    
                for(j=0;j<=B-1;j++)
                {                          
                        p=mypow(2,m-L)*j;  //用自己定义的乘方函数,效率比库函数高很多.          
                        ws.real=cos_tab[p];
                        ws.imag=sin_tab[p];
                        for(i=j;i<=N-1;i=i+le)//遍历M级所有的碟形
                        {
                                ip=i+B;                 
                                //执行复数乘法
                                t.real=xin[ip].real*ws.real-xin[ip].imag*ws.imag;
                                t.imag=xin[ip].real*ws.imag+xin[ip].imag*ws.real;

                                xin[ip].real=xin.real-t.real;
                                xin[ip].imag=xin.imag-t.imag;
                                xin.real=xin.real+t.real;
                                xin.imag=xin.imag+t.imag;
                        }
                }
        }   
}




沙发
734774645| | 2015-12-29 22:24 | 只看该作者
使用STM32 的DSP库进行FFT变换
/*
*********************************************************************************************************
FileName:dsp_asm.h
*********************************************************************************************************
*/

#ifndef  __DSP_ASM_H__
#define  __DSP_ASM_H__
*********************************************************************************************************
*                                            FUNCTION PROTOTYPES
*********************************************************************************************************
*/

void    dsp_asm_test(void);
void  dsp_asm_init(void);

#endif                                                          /* End of module include.                               */
/*8888888888888888888888888888888888888888888888888888888888888888*/
/*8888888888888888888888888888888888888888888888888888888888888888*/
/*
* FileName:dsp_asm.c
* Author:Bobby.Chen
* Email:heroxx@163.com
* Date:2010-08-11
* Description:This file showes how to use the dsp library in mdk project.
*                  使用三角函数生成采样点,供FFT计算
*                  进行FFT测试时,按下面顺序调用函数即可:
*                   dsp_asm_init();
*                   dsp_asm_test();
*/
#include "stm32f10x.h"
#include "dsp_asm.h"
#include "stm32_dsp.h"
#include "table_fft.h"
#include <stdio.h>
#include <math.h>


/*
*********************************************************************************************************
*                                           LOCAL CONSTANTS
*********************************************************************************************************
*/
#define PI2 6.28318530717959
// Comment the lines that you don't want to use.
// 要模拟FFT,请注释掉其他的预定义
// 此处也可以全部注释掉,在MDK的工程属性->"C/C++"->"Preprocessor Symbols"-"Define:"中添加NPT_XXX项目
// 但是这样做法的缺点是每次修改XXX数据,都会导致MDK下次编译时会编译全部文件,速度太慢。
//#define NPT_64 64
#define NPT_256 256
//#define NPT_1024 1024

// N=64,Fs/N=50Hz,Max(Valid)=1600Hz
// 64点FFt,采样率3200Hz,频率分辨率50Hz,测量最大有效频率1600Hz
#ifdef NPT_64
#define NPT 64
#define Fs  3200
#endif

// N=256,Fs/N=25Hz,Max(Valid)=3200Hz
// 256点FFt,采样率6400Hz,频率分辨率25Hz,测量最大有效频率3200Hz
#ifdef NPT_256
#define NPT 256
#define Fs  6400
#endif

// N=1024,Fs/N=5Hz,Max(Valid)=2560Hz
// 1024点FFt,采样率5120Hz,频率分辨率5Hz,测量最大有效频率2560Hz
#ifdef NPT_1024
#define NPT 1024
#define Fs  5120
#endif

/*
*********************************************************************************************************
*                                       LOCAL GLOBAL VARIABLES
*********************************************************************************************************
*/
extern  uint16_t  TableFFT[];
long lBUFIN[NPT];         /* Complex input vector */
long lBUFOUT[NPT];        /* Complex output vector */
long lBUFMAG[NPT];/* Magnitude vector */
/*
*********************************************************************************************************
*                                      LOCAL FUNCTION PROTOTYPES
*********************************************************************************************************
*/
void dsp_asm_powerMag(void);

/*
*********************************************************************************************************
*   Initialize data tables for lBUFIN
* 模拟采样数据,采样数据中包含3种频率正弦波:50Hz,2500Hz,2550Hz
* lBUFIN数组中,每个单元数据高字(高16位)中存储采样数据的实部,低字(低16位)存储采样数据的虚部(总是为0)
*********************************************************************************************************
*/
void  dsp_asm_init()
{
  u16 i=0;
  float fx;
  for(i=0;i<NPT;i++)
  {
    fx  = 4000 * sin(PI2*i*50.0/Fs) + 4000 * sin(PI2*i*2500.0/Fs) + 4000*sin(PI2*i*2550.0/Fs);
    lBUFIN = ((s16)fx)<<16;
  }
}

/*
*********************************************************************************************************
*   Test FFT,calculate powermag
*   进行FFT变换,并计算各次谐波幅值
*********************************************************************************************************
*/
void  dsp_asm_test()
{
// 根据预定义选择合适的FFT函数
#ifdef NPT_64
  cr4_fft_64_stm32(lBUFOUT, lBUFIN, NPT);
#endif

#ifdef NPT_256
  cr4_fft_256_stm32(lBUFOUT, lBUFIN, NPT);
#endif

#ifdef NPT_1024
  cr4_fft_1024_stm32(lBUFOUT, lBUFIN, NPT);
#endif

  // 计算幅值
  dsp_asm_powerMag();

//  printf("No. Freq      Power\n");
//  for(i=0;i<NPT/2;i++)
//  {
//    printf("%4d,%4d,%10d,%10d,%10d\n",i,(u16)((float)i*Fs/NPT),lBUFMAG,(lBUFOUT>>16),(lBUFOUT&0xffff));
//  }
//  printf("*********END**********\r\n");
}
/*
*********************************************************************************************************
*   Calculate powermag
*   计算各次谐波幅值
*   先将lBUFOUT分解成实部(X)和虚部(Y),然后计算赋值(sqrt(X*X+Y*Y)
*********************************************************************************************************
*/
void dsp_asm_powerMag(void)
{
  s16 lX,lY;
  u32 i;
  for(i=0;i<NPT/2;i++)
  {
    lX  = (lBUFOUT << 16) >> 16;
    lY  = (lBUFOUT >> 16);
    {
    float X    = NPT * ((float)lX) /32768;
    float Y    = NPT * ((float)lY) /32768;
    float Mag = sqrt(X*X + Y*Y)/NPT;
    lBUFMAG    = (u32)(Mag * 65536);
    }
  }
}


// 笔者使用的是金牛开发板,CPU为STM32F107VC;JLink V8,MDK-ARM 4.10
// 注意FFT运算结果的对称性,也即256点的运算结果,只有前面128点的数据是有效可用的。

上图中,数组下标X对应的谐波频率为:N×Fs/64=N×3200/64=N*50Hz.

lBUFMAG[1] 对应 50Hz谐波幅值

上图中由于FFT分辨率50HZ,最大只能识别1600Hz谐波,导致结果中出现错误的数据。
// 256点FFT运算结果图(局部):


使用特权

评论回复
板凳
734774645| | 2015-12-29 22:25 | 只看该作者

上图中,数组下标X对应的谐波频率为:N×Fs/256=N×6400/256=N*25Hz.

lBUFMAG[2] 对应 2×25 =50Hz谐波幅值

lBUFMAG[100] 对应 100×25=2500Hz谐波幅值

lBUFMAG[102] 对应 102×25=2550Hz谐波幅值


// 1024点FFT运算结果图(局部):

上图中,数组下标X对应的谐波频率为:N×Fs/1024=N×5120/1024=N*5Hz.

lBUFMAG[10] 对应 10×5 =50Hz谐波幅值

lBUFMAG[500] 对应 500×5=2500Hz谐波幅值

lBUFMAG[510] 对应 510×5=2550Hz谐波幅值


该工程中模拟信号源为:4000 * sin(PI2*i*50.0/Fs) + 4000 * sin(PI2*i*2500.0/Fs) + 4000*sin(PI2*i*2550.0/Fs)

信号为1个50Hz、1个2500Hz、1个2550Hz的正弦波混合信号,幅值为均为4000。

点击这里下载示例工程


使用特权

评论回复
地板
892953881|  楼主 | 2015-12-29 22:38 | 只看该作者
734774645 发表于 2015-12-29 22:25
上图中,数组下标X对应的谐波频率为:N×Fs/256=N×6400/256=N*25Hz.lBUFMAG[2] 对应 2×25 =50Hz谐波幅值 ...

真是非常感谢了!!!写的这么详细,只是有点遗憾你放的源码包无法下载。

使用特权

评论回复
5
892953881|  楼主 | 2016-1-2 12:40 | 只看该作者
734774645 发表于 2015-12-29 22:25
上图中,数组下标X对应的谐波频率为:N×Fs/256=N×6400/256=N*25Hz.lBUFMAG[2] 对应 2×25 =50Hz谐波幅值 ...

3Q,在你的帮助下问题解决了

使用特权

评论回复
6
734774645| | 2016-1-8 20:39 | 只看该作者
很高兴帮到了楼主,楼主有好东西也赶紧发来论坛分享分享啊

使用特权

评论回复
7
mintspring| | 2016-1-8 22:35 | 只看该作者
用自己定义的乘方函数,效率比库函数高很多.这里如果采用移位计算,效率更高.      

使用特权

评论回复
8
网络孤客| | 2016-1-10 14:28 | 只看该作者
谢谢,一直没搞懂,慢慢研究。

使用特权

评论回复
9
米尔豪斯| | 2016-1-10 15:14 | 只看该作者
收藏先,stm32的DSP库还没用过呢

使用特权

评论回复
10
892953881|  楼主 | 2016-1-10 21:30 | 只看该作者
对于我自己贴的那个FFT的函数,函数的输入是复数,有人告诉我说,把电压值放在实数部位,虚数为0带进去计算。不知有没有人知道不是是这样的。

使用特权

评论回复
11
冷火智能| | 2017-1-11 16:08 | 只看该作者
892953881 发表于 2016-1-10 21:30
对于我自己贴的那个FFT的函数,函数的输入是复数,有人告诉我说,把电压值放在实数部位,虚数为0带进去计算 ...

楼主~最近写程序也是需要FFT,不知道那个复数的结构体是怎么定义的呢?

使用特权

评论回复
发新帖 我要提问
您需要登录后才可以回帖 登录 | 注册

本版积分规则

13

主题

42

帖子

2

粉丝