FPGA

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

2019-07-15 20:51发布

站内问答 / FPGA
2916 3 1474
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"

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

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

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

  62.         else if (n > FFT_N / 4 && n < FFT_N / 2) {
  63.                 n -= FFT_N / 4;
  64.                 a = SIN_TAB[FFT_N / 4 - n];
  65.         }

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

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

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

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

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

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

  11. #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"


  23. #define soc_cv_av
  24. #include "hwlib.h"
  25. #include "socal/socal.h"
  26. #include "socal/hps.h"

  28. #include <linux/types.h>
  29. #include <linux/spi/spidev.h>

  31. #define MY_CDEV_MAGIC 'k'
  32. #define MY_CDEV_CMD  _IO(MY_CDEV_MAGIC,0x1a)
  33. #define MY_CDEV_CMD_READ  _IOR(MY_CDEV_MAGIC,0x1b,int)
  34. #define DEVICE_MAJOR 205

  36. #define BUFFER_SIZE 1024*16*4
  37. int8_t ch1_data[WAVE_SIZE]; //通道1波形存储器
  38. int8_t ch2_data[WAVE_SIZE]; //通道2波形存储器
  39. #define HW_REGS_BASE (ALT_STM_OFST )
  40. #define HW_REGS_SPAN (0x04000000 )
  41. #define HW_REGS_MASK (HW_REGS_SPAN - 1 )
  42. char *fb;


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

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


  52. float Amp[20] = { };
  53. double freq[20] = { };

  55. static volatile unsigned long *pio_state_addr = NULL;
  56. static volatile unsigned long *pio_en_addr = NULL;
  57. static volatile unsigned long *pio_length_addr = NULL;
  58. static volatile unsigned long *pio_ambase_addr = NULL;
  59. static volatile unsigned long *pio_sfreq_addr = NULL;
  60. static volatile unsigned long *pio_trig_addr = NULL;
  61. static volatile unsigned long *touch_addr = NULL;
  62. static volatile unsigned long *ad9226_measure_addr = NULL;
  63. static volatile unsigned long *pio_led_addr = NULL;

  65. int FPGA_init(void) {
  66.         int fd;
  67.         void *periph_virtual_base;
  68.         if ((fd = open("/dev/mem", ( O_RDWR | O_SYNC))) == -1) {
  69.                 printf("ERROR: could not open "/dev/mem"... ");
  70.                 return (1);
  71.         }
  72.         periph_virtual_base = mmap( NULL, HW_REGS_SPAN, ( PROT_READ | PROT_WRITE),
  73.         MAP_SHARED, fd, HW_REGS_BASE);
  74.         if (periph_virtual_base == MAP_FAILED) {
  75.                 printf("ERROR: mmap() failed... ");
  76.                 close(fd);
  77.                 return (1);
  78.         }
  79.         pio_state_addr = periph_virtual_base
  80.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_STATE_BASE)
  81.                                         & (unsigned long) ( HW_REGS_MASK));
  82.         pio_en_addr = periph_virtual_base
  83.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_EN_BASE)
  84.                                         & (unsigned long) ( HW_REGS_MASK));
  85.         pio_length_addr = periph_virtual_base
  86.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_LENGTH_BASE)
  87.                                         & (unsigned long) ( HW_REGS_MASK));
  88.         pio_ambase_addr = periph_virtual_base
  89.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_AMBASE_BASE)
  90.                                         & (unsigned long) ( HW_REGS_MASK));
  91.         pio_sfreq_addr = periph_virtual_base
  92.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_SFREQ_SET_BASE)
  93.                                         & (unsigned long) ( HW_REGS_MASK));
  94.         pio_trig_addr = periph_virtual_base
  95.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_TRIG_SET_BASE)
  96.                                         & (unsigned long) ( HW_REGS_MASK));
  97.         touch_addr = periph_virtual_base
  98.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + XPT2046_BASE)
  99.                                         & (unsigned long) ( HW_REGS_MASK));
  100.         ad9226_measure_addr = periph_virtual_base
  101.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + AD9226_MEASURE_BASE)
  102.                                         & (unsigned long) ( HW_REGS_MASK));
  103.         pio_led_addr = periph_virtual_base
  104.                         + ((unsigned long) ( ALT_LWFPGASLVS_OFST + LED_PIO_BASE)
  105.                                         & (unsigned long) ( HW_REGS_MASK));

  107.         printf("pio_led_addr is %x ", pio_led_addr);
  108.         return 0;
  109. }

  111. int init_fb(long int *screensize) {
  112.         int fbfd = 0;
  113.         struct fb_var_screeninfo vinfo;
  114.         struct fb_fix_screeninfo finfo;

  116.         fbfd = open("/dev/fb0", O_RDWR);
  117.         // Open the file for reading and writing
  118.         if (!fbfd) {
  119.                 printf("Error: cannot open framebuffer device. ");
  120.                 exit(0);
  121.         }
  122.         printf("The framebuffer device was opened successfully. ");

  124. // Get fixed screen information
  125.         if (ioctl(fbfd, FBIOGET_FSCREENINFO, &finfo)) {
  126.                 printf("Error reading fixed information. ");
  127.                 exit(0);
  128.         }

  130. // Get variable screen information
  131.         if (ioctl(fbfd, FBIOGET_VSCREENINFO, &vinfo)) {
  132.                 printf("Error reading variable information. ");
  133.                 exit(0);
  134.         }
  135.         printf("%dx%d, %dbpp ", vinfo.xres, vinfo.yres, vinfo.bits_per_pixel);
  136.         // Figure out the size of the screen in bytes
  137.         *screensize = vinfo.xres * vinfo.yres * vinfo.bits_per_pixel / 8;
  138.         // Map the device to memory
  139.         fb = (char *) mmap(0, *screensize, PROT_READ | PROT_WRITE, MAP_SHARED, fbfd,
  140.                         0);

  142.         if ((int) fb == -1) {
  143.                 printf("Error: failed to map framebuffer device to memory. ");
  144.                 exit(0);
  145.         }

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

  149.         return fbfd;
  150. }
  151. void osc_init() {
  152.         LCD_Dir(LCD_HORIZONTAL);
  153.         draw_wave_area();
  154.         draw_main_ui();

  156.         CH1_Y_Offset = 64;
  157.         CH2_Y_Offset = -64;
  158.         CH1_X_Offset = 20;
  159.         CH2_X_Offset = 0;

  161.         CH1_TIME_factor = 1;
  162.         CH2_TIME_factor = CH1_TIME_factor;
  163. }

  165. int sample_init(void) {
  166.         int fd;
  167.         unsigned long dma_base;
  168.         fd = open("/dev/de1_soc_demo", O_RDWR);
  169.         if (fd == 1) {
  170.                 perror("open error ");
  171.                 exit(-1);
  172.         }

  174.         ioctl(fd, MY_CDEV_CMD, &dma_base);
  175.         printf("dma_base is %x ", dma_base);

  177.         *pio_ambase_addr = dma_base;
  178.         *pio_length_addr = 16384;
  179.         *pio_trig_addr = 2047 | 0x8000;
  180.         *pio_sfreq_addr = 0xf;

  182.         return fd;
  183. }


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

  188.         int fd;
  189.         int fbfd = 0;

  191.         int32_t read_buffer[NPT] = { 0 };
  192.         int32_t retval;
  193.         long int screensize = 0;

  195.         int32_t dtmp[NPT] = { };

  197.         int32_t AppU, AppD;
  198.         float App;

  200.         int32_t i, j;
  201.         unsigned int tmp = 0;

  203.         float max_ch0 = 0;
  204.         int32_t max_num_ch0 = 0;


  207.         int32_t t;
  208.         float sample_rate;
  209.         sample_rate = 2048;
  210.         create_sin_tab(&my_sin_tab, NPT);


  213.         fpga_init();
  214.         fbfd = init_fb(&screensize);
  215.         osc_init();
  216.         fd = sample_init();

  218.         *pio_en_addr = 2;
  219.         *pio_en_addr = 0;

  221.         while (1) {
  222.                 *pio_sfreq_addr = 50000000/sample_rate - 1;
  223.                 *pio_en_addr = 1;
  224.                 *pio_en_addr = 0;
  225.                 while (!(*pio_state_addr))
  226.                         ;
  227.                 retval = read(fd, read_buffer, NPT);

  229.                 for(i = 0;i < NPT; i++ ){
  230.                                 tmp = read_buffer[i]&0x0000ffff;                        //chanel 0
  231.                                 if( (tmp&0x800) != 0){                                                //negative
  232.                                         tmp = (tmp | 0xfffff000);
  233.                                 }
  234.                                 dtmp[i] = (int)tmp;
  235.                         }

  237.                 AppU = -2048;
  238.                 AppD = 2047;
  239.                 for (i = 0; i < NPT; i++) {
  240.                         if (dtmp[i] > AppU)
  241.                                 AppU = dtmp[i];
  242.                         if (dtmp[i] < AppD)
  243.                                 AppD = dtmp[i];
  244.                 }

  246.                 for (t = 0; t < NPT; t++) {
  247.                         s[t].real = 0;
  248.                         s[t].imag = 0;
  249.                         ADC_ConvertedValueLocal = 0;
  250.                         i = 0;
  251.                         ADC_ConvertedValueLocal = (float) dtmp[t] * (5000.0 / 2048);
  252.                         s[t].real = ADC_ConvertedValueLocal;
  253.                 }
  254.                 FFT(s);

  256.                 for (i = 0; i < NPT; i++) { /* 求变换后结果的模值，存入复数的实部部分 */
  257.                         s[i].real = sqrt(s[i].real * s[i].real + s[i].imag * s[i].imag)
  258.                                         / (i == 0 ? NPT : NPT / 2);
  259.                 }
  260.                 max_ch0 = 0;
  261.                 max_num_ch0 = 0;
  262.                 for (i = 1; i < NPT/2; i++) {
  263.                         if (s[i].real > max_ch0) {
  264.                                 max_ch0 = s[i].real;
  265.                                 max_num_ch0 = i;
  266.                         }
  267.                 }

  269.                 freq[0] = sample_rate * 1.0 * max_num_ch0 / NPT;
  270.                 Amp[0] = max_ch0;

  272. //                for (i = 1; i < 20; i++) {
  273. //                        freq[i] = freq[0] * (i + 1);
  274. //                        Amp[i] = s[(i + 1) * max_num_ch0].real;
  275. //
  276. //                        printf("Amp[%d]:%.1fmA ",i, Amp[i]);
  277. //                        printf("freq[%d]:%.1fHz ",i, freq[i]);
  278. //                }
  279.                 App = 5000.0 * (AppU - AppD) / 2048;

  281.                 printf("App:%.1fmA ", App);
  282.                 printf("Amp:%.1fmA ", Amp[0]*4);
  283.                 printf("Fre:%.1fHz ", freq[0]);

  285.                 j++;
  286.                 printf("%d", j);


  289.                 if (retval == -1) {
  290.                         perror("read error ");
  291.                         exit(-1);
  292.                 }

  294.                 usleep(20000);

  296.         }
  297.         printf("The app is closed. ");
  298.         munmap(fb, screensize);
  299.         close(fbfd);
  300.         close(fd);
  301.         return 0;

  303. }
复制代码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。




友情提示: 此问题已得到解决，问题已经关闭,关闭后问题禁止继续编辑，回答。
该问题目前已经被作者或者管理员关闭, 无法添加新回复
3条回答
李春明
1楼-- · 2019-07-16 02:36
学习一下
liujinyi016
2楼-- · 2019-07-16 07:57
你的FFT有没有要求字节对齐呢
小梅哥
3楼-- · 2019-07-16 11:05
.问题找到了。main.c中227行read函数用错了,这里应该是读点数的4倍的字节数据，之前只读了点数个字节的数据，导致分析时后面3/4个周期的数据全为0，所以导致计算得到的幅度为理论值的1/4

一周热门 更多>