【2018电赛A题】关于Linux下做FFT运算的一疑惑

2019-07-15 20:51发布

2018年ti杯多省联赛落下帷幕。看了赛题,觉得还是蛮有意思的,第一个题目,是进行电流幅度的计算,需要用到FFT运算。我这里是在Intel Cyclone V SoC上进行了一次实验,使用12位有符号的高速ADC(AD9226),通过控制采样设置采样率为2048Hz,采样8192个点,然后执行FFT运算。在网上找了一个成熟的FFT运算代码,加入工程中编译,结果发现计算的幅度值永远只有实际幅度值的1/4,而该代码在NIOS II CPU上运行却是正确的。一直没找到问题所在。特发代码出来,看有没有人能指点一二。
FFT.c
  1. /*********************************************************************
  2. 快速福利叶变换C程序包
  3. 函数简介:此程序包是通用的快速傅里叶变换C语言函数,移植性强,以下部分不依
  4. 赖硬件。此程序包采用联合体的形式表示一个复数,输入为自然顺序的复
  5. 数(输入实数是可令复数虚部为0),输出为经过FFT变换的自然顺序的
  6. 复数.此程序包可在初始化时调用create_sin_tab()函数创建正弦函数表,
  7. 以后的可采用查表法计算耗时较多的sin和cos运算,加快可计算速度.与
  8. Ver1.1版相比较,Ver1.2版在创建正弦表时只建立了1/4个正弦波的采样值,
  9. 相比之下节省了FFT_N/4个存储空间
  10. 使用说明:使用此函数只需更改宏定义FFT_N的值即可实现点数的改变,FFT_N的
  11. 应该为2的N次方,不满足此条件时应在后面补0。若使用查表法计算sin值和
  12. cos值,应在调用FFT函数前调用create_sin_tab()函数创建正弦表
  13. 函数调用:FFT(s);
  14. 版    本:Ver1.2
  15. 参考文献:
  16. **********************************************************************/
  17. #include <math.h>
  18. #include "fft.h"

  19. float *SIN_TAB; //定义正弦表的存放空间
  20. int FFT_N = 8192; //定义采样点大小
  21. /*******************************************************************
  22. 函数原型:struct compx EE(struct compx b1,struct compx b2)
  23. 函数功能:对两个复数进行乘法运算
  24. 输入参数:两个以联合体定义的复数a,b
  25. 输出参数:a和b的乘积,以联合体的形式输出
  26. *******************************************************************/
  27. struct compx EE(struct compx a, struct compx b) {
  28.         struct compx c;
  29.         c.real = a.real * b.real - a.imag * b.imag;
  30.         c.imag = a.real * b.imag + a.imag * b.real;
  31.         return (c);
  32. }

  33. /******************************************************************
  34. 函数原型:void create_sin_tab(float *sin_t,int PointNum)
  35. 函数功能:创建一个正弦采样表,采样点数与福利叶变换点数相同
  36. 输入参数:*sin_t存放正弦表的数组指针,PointNum采样点数
  37. 输出参数:无
  38. ******************************************************************/
  39. void create_sin_tab(float *sin_t, int PointNum) {
  40.         int i;
  41.         SIN_TAB = sin_t;
  42.         FFT_N = PointNum;
  43.         for (i = 0; i <= FFT_N / 4; i++)
  44.                 SIN_TAB[i] = sin(2 * PI * i / FFT_N);
  45. }
  46. /******************************************************************
  47. 函数原型:void sin_tab(float pi)
  48. 函数功能:采用查表的方法计算一个数的正弦值
  49. 输入参数:pi 所要计算正弦值弧度值,范围0--2*PI,不满足时需要转换
  50. 输出参数:输入值pi的正弦值
  51. ******************************************************************/
  52. float sin_tab(float pi) {
  53.         int n = 0;
  54.         float a = 0;
  55.         n = (int) (pi * FFT_N / 2 / PI);

  56.         if (n >= 0 && n <= FFT_N / 4)
  57.                 a = SIN_TAB[n];

  58.         else if (n > FFT_N / 4 && n < FFT_N / 2) {
  59.                 n -= FFT_N / 4;
  60.                 a = SIN_TAB[FFT_N / 4 - n];
  61.         }

  62.         else if (n >= FFT_N / 2 && n < 3 * FFT_N / 4) {
  63.                 n -= FFT_N / 2;
  64.                 a = -SIN_TAB[n];
  65.         } else if (n >= 3 * FFT_N / 4 && n < 3 * FFT_N) {
  66.                 n = FFT_N - n;
  67.                 a = -SIN_TAB[n];
  68.         }
  69.         return a;
  70. }
  71. /******************************************************************
  72. 函数原型:void cos_tab(float pi)
  73. 函数功能:采用查表的方法计算一个数的余弦值
  74. 输入参数:pi 所要计算余弦值弧度值,范围0--2*PI,不满足时需要转换
  75. 输出参数:输入值pi的余弦值
  76. ******************************************************************/
  77. float cos_tab(float pi) {
  78.         float a, pi2;
  79.         pi2 = pi + PI / 2;
  80.         if (pi2 > 2 * PI)
  81.                 pi2 -= 2 * PI;
  82.         a = sin_tab(pi2);
  83.         return a;
  84. }
  85. /*****************************************************************
  86. 函数原型:void FFT(struct compx *xin)
  87. 函数功能:对输入的复数组进行快速傅里叶变换(FFT)
  88. 输入参数:*xin复数结构体组的首地址指针,struct型
  89. 输出参数:无
  90. *****************************************************************/
  91. void FFT(struct compx *xin) {
  92.         register int nv2, nm1;
  93.         struct compx u, w, t;
  94.         int f, m, i, k, l, j = 0;
  95.         nv2 = FFT_N / 2; //变址运算,即把自然顺序变成倒位序,采用雷德算法
  96.         nm1 = FFT_N - 1;

  97.         for (i = 0; i < nm1; i++) {
  98.                 if (i < j) { //如果i<j,即进行变址
  99.                         t = xin[j]; //
  100.                         xin[j] = xin[i]; //
  101.                         xin[i] = t;
  102.                 } //
  103.                 k = nv2; //求j的下一个倒位序
  104.                 while (k <= j) { //如果k<=j,表示j的最高位为1
  105.                         j = j - k; //把最高位变成0
  106.                         k = k / 2;
  107.                 } //k/2,比较次高位,依次类推,逐个比较,直到某个位为0
  108.                 j = j + k;
  109.         } //把0改为1
  110.         {
  111.                 int le, lei, ip; //FFT运算核,使用蝶形运算完成FFT运算
  112.                 f = FFT_N;
  113.                 for (l = 1; (f = f / 2) != 1; l++)
  114.                         ; //计算l的值,即计算蝶形级数
  115.                 for (m = 1; m <= l; m++) { // 控制蝶形结级数 //m表示第m级蝶形,l为蝶形级总数l=log(2)N
  116.                         le = 2 << (m - 1); //le蝶形结距离,即第m级蝶形的蝶形结相距le点
  117.                         lei = le / 2; //同一蝶形结中参加运算的两点的距离
  118.                         u.real = 1.0; //u为蝶形结运算系数,初始值为1
  119.                         u.imag = 0.0; //
  120.                         w.real = cos_tab(PI / lei); //w为系数商,即当前系数与前一个系数的商
  121.                         w.imag = -sin_tab(PI / lei); //

  122.                         for (j = 0; j <= lei - 1; j++) { //控制计算不同种蝶形结,即计算系数不同的蝶形结
  123.                                 for (i = j; i <= FFT_N - 1; i = i + le) { //控制同一蝶形结运算,即计算系数相同蝶形结
  124.                                         ip = i + lei; //i,ip分别表示参加蝶形运算的两个节点
  125.                                         t = EE(xin[ip], u); //蝶形运算,详见公式
  126.                                         xin[ip].real = xin[i].real - t.real; //
  127.                                         xin[ip].imag = xin[i].imag - t.imag; //
  128.                                         xin[i].real = xin[i].real + t.real; //
  129.                                         xin[i].imag = xin[i].imag + t.imag;
  130.                                 } //
  131.                                 u = EE(u, w);
  132.                         } //改变系数,进行下一个蝶形运算
  133.                 }
  134.         }
  135. }

复制代码FFT.h
  1. #ifndef FFT_H
  2. #define FFT_H
  3.                                                  //定义福利叶变换的点数
  4. #define PI 3.1415926535897932384626433832795028841971               //定义圆周率值

  5. struct compx {double real,imag;};                                    //定义一个复数结构

  6. extern void create_sin_tab(float *sin_t,int PointNum);
  7. extern void FFT(struct compx *xin);

  8. #endif // FFT_H
复制代码
main.c
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <unistd.h>
  4. #include <fcntl.h>
  5. #include <sys/mman.h>
  6. #include <sys/types.h>
  7. #include <dirent.h>
  8. #include <inttypes.h>
  9. #include <sys/time.h>
  10. #include <linux/fb.h>
  11. #include <stdbool.h>
  12. #include "stdint.h"
  13. #include "TFT_API.h"
  14. #include "data.h"
  15. #include "DSO.h"
  16. #include "test_data.h"
  17. #include <sys/ioctl.h>
  18. #include "math.h"
  19. #include "fft.h"
  20. #include "hps_0.h"


  21. #define soc_cv_av
  22. #include "hwlib.h"
  23. #include "socal/socal.h"
  24. #include "socal/hps.h"

  25. #include <linux/types.h>
  26. #include <linux/spi/spidev.h>

  27. #define MY_CDEV_MAGIC 'k'
  28. #define MY_CDEV_CMD  _IO(MY_CDEV_MAGIC,0x1a)
  29. #define MY_CDEV_CMD_READ  _IOR(MY_CDEV_MAGIC,0x1b,int)
  30. #define DEVICE_MAJOR 205

  31. #define BUFFER_SIZE 1024*16*4
  32. int8_t ch1_data[WAVE_SIZE]; //通道1波形存储器
  33. int8_t ch2_data[WAVE_SIZE]; //通道2波形存储器
  34. #define HW_REGS_BASE (ALT_STM_OFST )
  35. #define HW_REGS_SPAN (0x04000000 )
  36. #define HW_REGS_MASK (HW_REGS_SPAN - 1 )
  37. char *fb;


  38. #define                NPT          8192        /*FFT采样点数,该点数应该为2的N次方,不满足此条件时应在后面补0*/

  39. float my_sin_tab[NPT / 4 + 1]; /*定义正弦表的存放空间*/
  40. struct compx s[NPT]; /*FFT输入和输出:从S[0]开始存放,根据大小自己定义*/
  41. float ADC_ConvertedValueLocal; /*自己定义的局部变量,用于保存转换计算后的电压值         */


  42. float Amp[20] = { };
  43. double freq[20] = { };

  44. static volatile unsigned long *pio_state_addr = NULL;
  45. static volatile unsigned long *pio_en_addr = NULL;
  46. static volatile unsigned long *pio_length_addr = NULL;
  47. static volatile unsigned long *pio_ambase_addr = NULL;
  48. static volatile unsigned long *pio_sfreq_addr = NULL;
  49. static volatile unsigned long *pio_trig_addr = NULL;
  50. static volatile unsigned long *touch_addr = NULL;
  51. static volatile unsigned long *ad9226_measure_addr = NULL;
  52. static volatile unsigned long *pio_led_addr = NULL;

  53. int FPGA_init(void) {
  54.         int fd;
  55.         void *periph_virtual_base;
  56.         if ((fd = open("/dev/mem", ( O_RDWR | O_SYNC))) == -1) {
  57.                 printf("ERROR: could not open "/dev/mem"... ");
  58.                 return (1);
  59.         }
  60.         periph_virtual_base = mmap( NULL, HW_REGS_SPAN, ( PROT_READ | PROT_WRITE),
  61.         MAP_SHARED, fd, HW_REGS_BASE);
  62.         if (periph_virtual_base == MAP_FAILED) {
  63.                 printf("ERROR: mmap() failed... ");
  64.                 close(fd);
  65.                 return (1);
  66.         }
  67.         pio_state_addr = periph_virtual_base
  68.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_STATE_BASE)
  69.                                         & (unsigned long) ( HW_REGS_MASK));
  70.         pio_en_addr = periph_virtual_base
  71.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_EN_BASE)
  72.                                         & (unsigned long) ( HW_REGS_MASK));
  73.         pio_length_addr = periph_virtual_base
  74.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_LENGTH_BASE)
  75.                                         & (unsigned long) ( HW_REGS_MASK));
  76.         pio_ambase_addr = periph_virtual_base
  77.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_AMBASE_BASE)
  78.                                         & (unsigned long) ( HW_REGS_MASK));
  79.         pio_sfreq_addr = periph_virtual_base
  80.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_SFREQ_SET_BASE)
  81.                                         & (unsigned long) ( HW_REGS_MASK));
  82.         pio_trig_addr = periph_virtual_base
  83.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_TRIG_SET_BASE)
  84.                                         & (unsigned long) ( HW_REGS_MASK));
  85.         touch_addr = periph_virtual_base
  86.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + XPT2046_BASE)
  87.                                         & (unsigned long) ( HW_REGS_MASK));
  88.         ad9226_measure_addr = periph_virtual_base
  89.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + AD9226_MEASURE_BASE)
  90.                                         & (unsigned long) ( HW_REGS_MASK));
  91.         pio_led_addr = periph_virtual_base
  92.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + LED_PIO_BASE)
  93.                                         & (unsigned long) ( HW_REGS_MASK));

  94.         printf("pio_led_addr is %x ", pio_led_addr);
  95.         return 0;
  96. }

  97. int init_fb(long int *screensize) {
  98.         int fbfd = 0;
  99.         struct fb_var_screeninfo vinfo;
  100.         struct fb_fix_screeninfo finfo;

  101.         fbfd = open("/dev/fb0", O_RDWR);
  102.         // Open the file for reading and writing
  103.         if (!fbfd) {
  104.                 printf("Error: cannot open framebuffer device. ");
  105.                 exit(0);
  106.         }
  107.         printf("The framebuffer device was opened successfully. ");

  108. // Get fixed screen information
  109.         if (ioctl(fbfd, FBIOGET_FSCREENINFO, &finfo)) {
  110.                 printf("Error reading fixed information. ");
  111.                 exit(0);
  112.         }

  113. // Get variable screen information
  114.         if (ioctl(fbfd, FBIOGET_VSCREENINFO, &vinfo)) {
  115.                 printf("Error reading variable information. ");
  116.                 exit(0);
  117.         }
  118.         printf("%dx%d, %dbpp ", vinfo.xres, vinfo.yres, vinfo.bits_per_pixel);
  119.         // Figure out the size of the screen in bytes
  120.         *screensize = vinfo.xres * vinfo.yres * vinfo.bits_per_pixel / 8;
  121.         // Map the device to memory
  122.         fb = (char *) mmap(0, *screensize, PROT_READ | PROT_WRITE, MAP_SHARED, fbfd,
  123.                         0);

  124.         if ((int) fb == -1) {
  125.                 printf("Error: failed to map framebuffer device to memory. ");
  126.                 exit(0);
  127.         }

  128.         printf("The framebuffer device was mapped to memory successfully. ");

  129.         return fbfd;
  130. }
  131. void osc_init() {
  132.         LCD_Dir(LCD_HORIZONTAL);
  133.         draw_wave_area();
  134.         draw_main_ui();

  135.         CH1_Y_Offset = 64;
  136.         CH2_Y_Offset = -64;
  137.         CH1_X_Offset = 20;
  138.         CH2_X_Offset = 0;

  139.         CH1_TIME_factor = 1;
  140.         CH2_TIME_factor = CH1_TIME_factor;
  141. }

  142. int sample_init(void) {
  143.         int fd;
  144.         unsigned long dma_base;
  145.         fd = open("/dev/de1_soc_demo", O_RDWR);
  146.         if (fd == 1) {
  147.                 perror("open error ");
  148.                 exit(-1);
  149.         }

  150.         ioctl(fd, MY_CDEV_CMD, &dma_base);
  151.         printf("dma_base is %x ", dma_base);

  152.         *pio_ambase_addr = dma_base;
  153.         *pio_length_addr = 16384;
  154.         *pio_trig_addr = 2047 | 0x8000;
  155.         *pio_sfreq_addr = 0xf;

  156.         return fd;
  157. }


  158. int main(int argc, char ** argv) {

  159.         int fd;
  160.         int fbfd = 0;

  161.         int32_t read_buffer[NPT] = { 0 };
  162.         int32_t retval;
  163.         long int screensize = 0;

  164.         int32_t dtmp[NPT] = { };

  165.         int32_t AppU, AppD;
  166.         float App;

  167.         int32_t i, j;
  168.         unsigned int tmp = 0;

  169.         float max_ch0 = 0;
  170.         int32_t max_num_ch0 = 0;


  171.         int32_t t;
  172.         float sample_rate;
  173.         sample_rate = 2048;
  174.         create_sin_tab(&my_sin_tab, NPT);


  175.         fpga_init();
  176.         fbfd = init_fb(&screensize);
  177.         osc_init();
  178.         fd = sample_init();

  179.         *pio_en_addr = 2;
  180.         *pio_en_addr = 0;

  181.         while (1) {
  182.                 *pio_sfreq_addr = 50000000/sample_rate - 1;
  183.                 *pio_en_addr = 1;
  184.                 *pio_en_addr = 0;
  185.                 while (!(*pio_state_addr))
  186.                         ;
  187.                 retval = read(fd, read_buffer, NPT);

  188.                 for(i = 0;i < NPT; i++ ){
  189.                                 tmp = read_buffer[i]&0x0000ffff;                        //chanel 0
  190.                                 if( (tmp&0x800) != 0){                                                //negative
  191.                                         tmp = (tmp | 0xfffff000);
  192.                                 }
  193.                                 dtmp[i] = (int)tmp;
  194.                         }

  195.                 AppU = -2048;
  196.                 AppD = 2047;
  197.                 for (i = 0; i < NPT; i++) {
  198.                         if (dtmp[i] > AppU)
  199.                                 AppU = dtmp[i];
  200.                         if (dtmp[i] < AppD)
  201.                                 AppD = dtmp[i];
  202.                 }

  203.                 for (t = 0; t < NPT; t++) {
  204.                         s[t].real = 0;
  205.                         s[t].imag = 0;
  206.                         ADC_ConvertedValueLocal = 0;
  207.                         i = 0;
  208.                         ADC_ConvertedValueLocal = (float) dtmp[t] * (5000.0 / 2048);
  209.                         s[t].real = ADC_ConvertedValueLocal;
  210.                 }
  211.                 FFT(s);

  212.                 for (i = 0; i < NPT; i++) { /* 求变换后结果的模值,存入复数的实部部分 */
  213.                         s[i].real = sqrt(s[i].real * s[i].real + s[i].imag * s[i].imag)
  214.                                         / (i == 0 ? NPT : NPT / 2);
  215.                 }
  216.                 max_ch0 = 0;
  217.                 max_num_ch0 = 0;
  218.                 for (i = 1; i < NPT/2; i++) {
  219.                         if (s[i].real > max_ch0) {
  220.                                 max_ch0 = s[i].real;
  221.                                 max_num_ch0 = i;
  222.                         }
  223.                 }

  224.                 freq[0] = sample_rate * 1.0 * max_num_ch0 / NPT;
  225.                 Amp[0] = max_ch0;

  226. //                for (i = 1; i < 20; i++) {
  227. //                        freq[i] = freq[0] * (i + 1);
  228. //                        Amp[i] = s[(i + 1) * max_num_ch0].real;
  229. //
  230. //                        printf("Amp[%d]:%.1fmA ",i, Amp[i]);
  231. //                        printf("freq[%d]:%.1fHz ",i, freq[i]);
  232. //                }
  233.                 App = 5000.0 * (AppU - AppD) / 2048;

  234.                 printf("App:%.1fmA ", App);
  235.                 printf("Amp:%.1fmA ", Amp[0]*4);
  236.                 printf("Fre:%.1fHz ", freq[0]);

  237.                 j++;
  238.                 printf("%d", j);


  239.                 if (retval == -1) {
  240.                         perror("read error ");
  241.                         exit(-1);
  242.                 }

  243.                 usleep(20000);

  244.         }
  245.         printf("The app is closed. ");
  246.         munmap(fb, screensize);
  247.         close(fbfd);
  248.         close(fd);
  249.         return 0;

  250. }
复制代码main文件里有些许冗余代码没有删掉,请直接看main函数中相关内容即可。
以下为完整的DS-5工程文件源码。 ADC_FFT.rar (83.11 KB, 下载次数: 74)
在运行结果中,如果对FFT的运算结果乘以4,得到的结果就是正确的幅值。如下图是乘以4之后的结果,一个是10Hz频率时候,一个是100Hz频率时候: 1.png 4.png
再次说明下参数:2048Hz采样率,8192个点FFT,ADC是12位有符号ADC。使用的DMA完成数据的采集,DMA采集了2个通道共24位的数据,然后高16位存储通道1的数据,低16位存储通道0的数据,需要注意的是。每个16位的数据中,低12位为数据,高12位始终为0。




友情提示: 此问题已得到解决,问题已经关闭,关闭后问题禁止继续编辑,回答。