` 2018年ti杯多省联赛落下帷幕。看了赛题,觉得还是蛮有意思的,第一个题目,是进行电流幅度的计算,需要用到FFT运算。我这里是在Intel Cyclone V SoC上进行了一次实验,使用12位有符号的高速ADC(AD9226),通过控制采样设置采样率为2048Hz,采样8192个点,然后执行FFT运算。在网上找了一个成熟的FFT运算代码,加入工程中编译,结果发现计算的幅度值永远只有实际幅度值的1/4,而该代码在NIOS II CPU上运行却是正确的。一直没找到问题所在。特发代码出来,看有没有人能指点一二。
FFT.c
- /*********************************************************************
- 快速福利叶变换C程序包
- 函数简介:此程序包是通用的快速傅里叶变换C语言函数,移植性强,以下部分不依
- 赖硬件。此程序包采用联合体的形式表示一个复数,输入为自然顺序的复
- 数(输入实数是可令复数虚部为0),输出为经过FFT变换的自然顺序的
- 复数.此程序包可在初始化时调用create_sin_tab()函数创建正弦函数表,
- 以后的可采用查表法计算耗时较多的sin和cos运算,加快可计算速度.与
- Ver1.1版相比较,Ver1.2版在创建正弦表时只建立了1/4个正弦波的采样值,
- 相比之下节省了FFT_N/4个存储空间
- 使用说明:使用此函数只需更改宏定义FFT_N的值即可实现点数的改变,FFT_N的
- 应该为2的N次方,不满足此条件时应在后面补0。若使用查表法计算sin值和
- cos值,应在调用FFT函数前调用create_sin_tab()函数创建正弦表
- 函数调用:FFT(s);
- 版 本:Ver1.2
- 参考文献:
- **********************************************************************/
- #include
- #include "fft.h"
- float *SIN_TAB; //定义正弦表的存放空间
- int FFT_N = 8192; //定义采样点大小
- /*******************************************************************
- 函数原型:struct compx EE(struct compx b1,struct compx b2)
- 函数功能:对两个复数进行乘法运算
- 输入参数:两个以联合体定义的复数a,b
- 输出参数:a和b的乘积,以联合体的形式输出
- *******************************************************************/
- struct compx EE(struct compx a, struct compx b) {
- struct compx c;
- c.real = a.real * b.real - a.imag * b.imag;
- c.imag = a.real * b.imag + a.imag * b.real;
- return (c);
- }
- /******************************************************************
- 函数原型:void create_sin_tab(float *sin_t,int PointNum)
- 函数功能:创建一个正弦采样表,采样点数与福利叶变换点数相同
- 输入参数:*sin_t存放正弦表的数组指针,PointNum采样点数
- 输出参数:无
- ******************************************************************/
- void create_sin_tab(float *sin_t, int PointNum) {
- int i;
- SIN_TAB = sin_t;
- FFT_N = PointNum;
- for (i = 0; i <= FFT_N / 4; i++)
- SIN_TAB[i] = sin(2 * PI * i / FFT_N);
- }
- /******************************************************************
- 函数原型:void sin_tab(float pi)
- 函数功能:采用查表的方法计算一个数的正弦值
- 输入参数:pi 所要计算正弦值弧度值,范围0--2*PI,不满足时需要转换
- 输出参数:输入值pi的正弦值
- ******************************************************************/
- float sin_tab(float pi) {
- int n = 0;
- float a = 0;
- n = (int) (pi * FFT_N / 2 / PI);
- if (n >= 0 && n <= FFT_N / 4)
- a = SIN_TAB[n];
- else if (n > FFT_N / 4 && n < FFT_N / 2) {
- n -= FFT_N / 4;
- a = SIN_TAB[FFT_N / 4 - n];
- }
- else if (n >= FFT_N / 2 && n < 3 * FFT_N / 4) {
- n -= FFT_N / 2;
- a = -SIN_TAB[n];
- } else if (n >= 3 * FFT_N / 4 && n < 3 * FFT_N) {
- n = FFT_N - n;
- a = -SIN_TAB[n];
- }
- return a;
- }
- /******************************************************************
- 函数原型:void cos_tab(float pi)
- 函数功能:采用查表的方法计算一个数的余弦值
- 输入参数:pi 所要计算余弦值弧度值,范围0--2*PI,不满足时需要转换
- 输出参数:输入值pi的余弦值
- ******************************************************************/
- float cos_tab(float pi) {
- float a, pi2;
- pi2 = pi + PI / 2;
- if (pi2 > 2 * PI)
- pi2 -= 2 * PI;
- a = sin_tab(pi2);
- return a;
- }
- /*****************************************************************
- 函数原型:void FFT(struct compx *xin)
- 函数功能:对输入的复数组进行快速傅里叶变换(FFT)
- 输入参数:*xin复数结构体组的首地址指针,struct型
- 输出参数:无
- *****************************************************************/
- void FFT(struct compx *xin) {
- register int nv2, nm1;
- struct compx u, w, t;
- int f, m, i, k, l, j = 0;
- nv2 = FFT_N / 2; //变址运算,即把自然顺序变成倒位序,采用雷德算法
- nm1 = FFT_N - 1;
- for (i = 0; i < nm1; i++) {
- if (i < j) { //如果i
- t = xin[j]; //
- xin[j] = xin[i]; //
- xin[i] = t;
- } //
- k = nv2; //求j的下一个倒位序
- while (k <= j) { //如果k<=j,表示j的最高位为1
- j = j - k; //把最高位变成0
- k = k / 2;
- } //k/2,比较次高位,依次类推,逐个比较,直到某个位为0
- j = j + k;
- } //把0改为1
- {
- int le, lei, ip; //FFT运算核,使用蝶形运算完成FFT运算
- f = FFT_N;
- for (l = 1; (f = f / 2) != 1; l++)
- ; //计算l的值,即计算蝶形级数
- for (m = 1; m <= l; m++) { // 控制蝶形结级数 //m表示第m级蝶形,l为蝶形级总数l=log(2)N
- le = 2 << (m - 1); //le蝶形结距离,即第m级蝶形的蝶形结相距le点
- lei = le / 2; //同一蝶形结中参加运算的两点的距离
- u.real = 1.0; //u为蝶形结运算系数,初始值为1
- u.imag = 0.0; //
- w.real = cos_tab(PI / lei); //w为系数商,即当前系数与前一个系数的商
- w.imag = -sin_tab(PI / lei); //
- for (j = 0; j <= lei - 1; j++) { //控制计算不同种蝶形结,即计算系数不同的蝶形结
- for (i = j; i <= FFT_N - 1; i = i + le) { //控制同一蝶形结运算,即计算系数相同蝶形结
- ip = i + lei; //i,ip分别表示参加蝶形运算的两个节点
- t = EE(xin[ip], u); //蝶形运算,详见公式
- xin[ip].real = xin[i].real - t.real; //
- xin[ip].imag = xin[i].imag - t.imag; //
- xin[i].real = xin[i].real + t.real; //
- xin[i].imag = xin[i].imag + t.imag;
- } //
- u = EE(u, w);
- } //改变系数,进行下一个蝶形运算
- }
- }
- }
复制代码
FFT.h
- #ifndef FFT_H
- #define FFT_H
- //定义福利叶变换的点数
- #define PI 3.1415926535897932384626433832795028841971 //定义圆周率值
- struct compx {double real,imag;}; //定义一个复数结构
- extern void create_sin_tab(float *sin_t,int PointNum);
- extern void FFT(struct compx *xin);
- #endif // FFT_H
复制代码
main.c
- #include
- #include
- #include
- #include
- #include
- #include
- #include
- #include
- #include
- #include
- #include
- #include "stdint.h"
- #include "TFT_API.h"
- #include "data.h"
- #include "DSO.h"
- #include "test_data.h"
- #include
- #include "math.h"
- #include "fft.h"
- #include "hps_0.h"
- #define soc_cv_av
- #include "hwlib.h"
- #include "socal/socal.h"
- #include "socal/hps.h"
- #include
- #include
- #define MY_CDEV_MAGIC 'k'
- #define MY_CDEV_CMD _IO(MY_CDEV_MAGIC,0x1a)
- #define MY_CDEV_CMD_READ _IOR(MY_CDEV_MAGIC,0x1b,int)
- #define DEVICE_MAJOR 205
- #define BUFFER_SIZE 1024*16*4
- int8_t ch1_data[WAVE_SIZE]; //通道1波形存储器
- int8_t ch2_data[WAVE_SIZE]; //通道2波形存储器
- #define HW_REGS_BASE (ALT_STM_OFST )
- #define HW_REGS_SPAN (0x04000000 )
- #define HW_REGS_MASK (HW_REGS_SPAN - 1 )
- char *fb;
- #define NPT 8192 /*FFT采样点数,该点数应该为2的N次方,不满足此条件时应在后面补0*/
- float my_sin_tab[NPT / 4 + 1]; /*定义正弦表的存放空间*/
- struct compx s[NPT]; /*FFT输入和输出:从S[0]开始存放,根据大小自己定义*/
- float ADC_ConvertedValueLocal; /*自己定义的局部变量,用于保存转换计算后的电压值 */
- float Amp[20] = { };
- double freq[20] = { };
- static volatile unsigned long *pio_state_addr = NULL;
- static volatile unsigned long *pio_en_addr = NULL;
- static volatile unsigned long *pio_length_addr = NULL;
- static volatile unsigned long *pio_ambase_addr = NULL;
- static volatile unsigned long *pio_sfreq_addr = NULL;
- static volatile unsigned long *pio_trig_addr = NULL;
- static volatile unsigned long *touch_addr = NULL;
- static volatile unsigned long *ad9226_measure_addr = NULL;
- static volatile unsigned long *pio_led_addr = NULL;
- int FPGA_init(void) {
- int fd;
- void *periph_virtual_base;
- if ((fd = open("/dev/mem", ( O_RDWR | O_SYNC))) == -1) {
- printf("ERROR: could not open "/dev/mem"...
- ");
- return (1);
- }
- periph_virtual_base = mmap( NULL, HW_REGS_SPAN, ( PROT_READ | PROT_WRITE),
- MAP_SHARED, fd, HW_REGS_BASE);
- if (periph_virtual_base == MAP_FAILED) {
- printf("ERROR: mmap() failed...
- ");
- close(fd);
- return (1);
- }
- pio_state_addr = periph_virtual_base
- + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_STATE_BASE)
- & (unsigned long) ( HW_REGS_MASK));
- pio_en_addr = periph_virtual_base
- + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_EN_BASE)
- & (unsigned long) ( HW_REGS_MASK));
- pio_length_addr = periph_virtual_base
- + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_LENGTH_BASE)
- & (unsigned long) ( HW_REGS_MASK));
- pio_ambase_addr = periph_virtual_base
- + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_AMBASE_BASE)
- & (unsigned long) ( HW_REGS_MASK));
- pio_sfreq_addr = periph_virtual_base
- + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_SFREQ_SET_BASE)
- & (unsigned long) ( HW_REGS_MASK));
- pio_trig_addr = periph_virtual_base
- + ((unsigned long) ( ALT_LWFPGASLVS_OFST + PIO_TRIG_SET_BASE)
- & (unsigned long) ( HW_REGS_MASK));
- touch_addr = periph_virtual_base
- + ((unsigned long) ( ALT_LWFPGASLVS_OFST + XPT2046_BASE)
- & (unsigned long) ( HW_REGS_MASK));
- ad9226_measure_addr = periph_virtual_base
- + ((unsigned long) ( ALT_LWFPGASLVS_OFST + AD9226_MEASURE_BASE)
- & (unsigned long) ( HW_REGS_MASK));
- pio_led_addr = periph_virtual_base
- + ((unsigned long) ( ALT_LWFPGASLVS_OFST + LED_PIO_BASE)
- & (unsigned long) ( HW_REGS_MASK));
- printf("pio_led_addr is %x
- ", pio_led_addr);
- return 0;
- }
- int init_fb(long int *screensize) {
- int fbfd = 0;
- struct fb_var_screeninfo vinfo;
- struct fb_fix_screeninfo finfo;
- fbfd = open("/dev/fb0", O_RDWR);
- // Open the file for reading and writing
- if (!fbfd) {
- printf("Error: cannot open framebuffer device.
- ");
- exit(0);
- }
- printf("The framebuffer device was opened successfully.
- ");
- // Get fixed screen information
- if (ioctl(fbfd, FBIOGET_FSCREENINFO, &finfo)) {
- printf("Error reading fixed information.
- ");
- exit(0);
- }
- // Get variable screen information
- if (ioctl(fbfd, FBIOGET_VSCREENINFO, &vinfo)) {
- printf("Error reading variable information.
- ");
- exit(0);
- }
- printf("%dx%d, %dbpp
- ", vinfo.xres, vinfo.yres, vinfo.bits_per_pixel);
- // Figure out the size of the screen in bytes
- *screensize = vinfo.xres * vinfo.yres * vinfo.bits_per_pixel / 8;
- // Map the device to memory
- fb = (char *) mmap(0, *screensize, PROT_READ | PROT_WRITE, MAP_SHARED, fbfd,
- 0);
- if ((int) fb == -1) {
- printf("Error: failed to map framebuffer device to memory.
- ");
- exit(0);
- }
- printf("The framebuffer device was mapped to memory successfully.
- ");
- return fbfd;
- }
- void osc_init() {
- LCD_Dir(LCD_HORIZONTAL);
- draw_wave_area();
- draw_main_ui();
- CH1_Y_Offset = 64;
- CH2_Y_Offset = -64;
- CH1_X_Offset = 20;
- CH2_X_Offset = 0;
- CH1_TIME_factor = 1;
- CH2_TIME_factor = CH1_TIME_factor;
- }
- int sample_init(void) {
- int fd;
- unsigned long dma_base;
- fd = open("/dev/de1_soc_demo", O_RDWR);
- if (fd == 1) {
- perror("open error
- ");
- exit(-1);
- }
- ioctl(fd, MY_CDEV_CMD, &dma_base);
- printf("dma_base is %x
- ", dma_base);
- *pio_ambase_addr = dma_base;
- *pio_length_addr = 16384;
- *pio_trig_addr = 2047 | 0x8000;
- *pio_sfreq_addr = 0xf;
- return fd;
- }
- int main(int argc, char ** argv) {
- int fd;
- int fbfd = 0;
- int32_t read_buffer[NPT] = { 0 };
- int32_t retval;
- long int screensize = 0;
- int32_t dtmp[NPT] = { };
- int32_t AppU, AppD;
- float App;
- int32_t i, j;
- unsigned int tmp = 0;
- float max_ch0 = 0;
- int32_t max_num_ch0 = 0;
- int32_t t;
- float sample_rate;
- sample_rate = 2048;
- create_sin_tab(&my_sin_tab, NPT);
- fpga_init();
- fbfd = init_fb(&screensize);
- osc_init();
- fd = sample_init();
- *pio_en_addr = 2;
- *pio_en_addr = 0;
- while (1) {
- *pio_sfreq_addr = 50000000/sample_rate - 1;
- *pio_en_addr = 1;
- *pio_en_addr = 0;
- while (!(*pio_state_addr))
- ;
- retval = read(fd, read_buffer, NPT);
- for(i = 0;i < NPT; i++ ){
- tmp = read_buffer[i]&0x0000ffff; //chanel 0
- if( (tmp&0x800) != 0){ //negative
- tmp = (tmp | 0xfffff000);
- }
- dtmp[i] = (int)tmp;
- }
- AppU = -2048;
- AppD = 2047;
- for (i = 0; i < NPT; i++) {
- if (dtmp[i] > AppU)
- AppU = dtmp[i];
- if (dtmp[i] < AppD)
- AppD = dtmp[i];
- }
- for (t = 0; t < NPT; t++) {
- s[t].real = 0;
- s[t].imag = 0;
- ADC_ConvertedValueLocal = 0;
- i = 0;
- ADC_ConvertedValueLocal = (float) dtmp[t] * (5000.0 / 2048);
- s[t].real = ADC_ConvertedValueLocal;
- }
- FFT(s);
- for (i = 0; i < NPT; i++) { /* 求变换后结果的模值,存入复数的实部部分 */
- s[i].real = sqrt(s[i].real * s[i].real + s[i].imag * s[i].imag)
- / (i == 0 ? NPT : NPT / 2);
- }
- max_ch0 = 0;
- max_num_ch0 = 0;
- for (i = 1; i < NPT/2; i++) {
- if (s[i].real > max_ch0) {
- max_ch0 = s[i].real;
- max_num_ch0 = i;
- }
- }
- freq[0] = sample_rate * 1.0 * max_num_ch0 / NPT;
- Amp[0] = max_ch0;
- // for (i = 1; i < 20; i++) {
- // freq[i] = freq[0] * (i + 1);
- // Amp[i] = s[(i + 1) * max_num_ch0].real;
- //
- // printf("Amp[%d]:%.1fmA ",i, Amp[i]);
- // printf("freq[%d]:%.1fHz
- ",i, freq[i]);
- // }
- App = 5000.0 * (AppU - AppD) / 2048;
- printf("App:%.1fmA ", App);
- printf("Amp:%.1fmA ", Amp[0]*4);
- printf("Fre:%.1fHz
- ", freq[0]);
- j++;
- printf("%d", j);
- if (retval == -1) {
- perror("read error
- ");
- exit(-1);
- }
- usleep(20000);
- }
- printf("The app is closed.
- ");
- munmap(fb, screensize);
- close(fbfd);
- close(fd);
- return 0;
- }
复制代码
main文件里有些许冗余代码没有删掉,请直接看main函数中相关内容即可。
以下为完整的DS-5工程文件源码。
在运行结果中,如果对FFT的运算结果乘以4,得到的结果就是正确的幅值。如下图是乘以4之后的结果,一个是10Hz频率时候,一个是100Hz频率时候:
再次说明下参数:2048Hz采样率,8192个点FFT,ADC是12位有符号ADC。使用的DMA完成数据的采集,DMA采集了2个通道共24位的数据,然后高16位存储通道1的数据,低16位存储通道0的数据,需要注意的是。每个16位的数据中,低12位为数据,高12位始终为0。
`
1
|
3个回答
|
|
|