可用于CH32V103R8T6的快速傅里叶FFT变换库
本帖最后由 catnull 于 2021-8-30 18:55 编辑一、引言
很幸运申请到了南京沁恒自主设计的RISC-V内核CH32V103R8T6。它有72Mhz的系统主频，SRAM有20KB, FLASH有64KB也是够用了。开发板很精致，也设置了Arduino的接口，方便使用一些Arduino接口板开发板。上面自带了一个WCH-LINK 下载仿真器，可以很方便的进行程序的调试。片上丰富的资源，可以展开DSP的实验，特别是离散傅里叶变换的实验。这种数字信号处理的应用相当广泛，例如一些市面上的调音器，还有心率识别，等等，都可以使用傅里叶变化对时域信号进行处理，之后变成频域信号。
二、思路简介
网络上已经有了STM32官方的快速傅里叶变换的库了，相关的应用也开展了不少。但是咱们的CH32V103行不行？直接移植过来，因为涉及到底层的运算优化，非我的能力所能达到。但是作为实验的题目，先解决零到壹的问题。从无到有来解决最好了。
后面找到一本数字信号的书，名字叫做《实用数字信号处理——从原理到应用》， Digital Signal Processing A Practical Guide for Engineers and Scientists , 是一个美国人写的， Steven W. Smith。由北航的两位老师翻译成中文。认真读了一番，因为之前数学上恶补了关于傅里叶变换的知识，啃下了物理系的《数学物理方法》中的复变函数论，所以读起来还算有所收获。
参考里面的BASIC程序，我移植为C语言。先在QT上调试了一番，后面移植到了开发板上，算是达到能用的水平。
但是特别提醒一点，需要打开MountainRiver山河编译器里面的一些编译选项，才能支持数学库里面的正弦余弦运算，以及浮点数的模拟运算的功能（RISC-IMAC是没有实现浮点指令集的，浮点运算需要模拟为整数进行）。
打钩用于一些浮点数的输入输出，方便调试。
下图中的m 编译时时 gcc会加上 -lm，会调用数学库。没有加的或无法计算正余弦函数。
3、快速傅里叶变换FFT库文件中的代码
具体的代码如下：
头文件： cr2_fft.h
/******************
*
* Fast Fourier Transformer
*
* Author : 张远东莆田seighbang@126.com
*
*
********************/
#ifndef CR2_FFT_H
#define CR2_FFT_H
#define PI 3.14159265
#define FFT_N1024
typedef struct complex
{
float real ;
float imag ;
} complex;
int fft(complex* x);
int ifft(complex* x );
#endif // CR2_FFT_H
库文件: cr2_fft.c#include "cr2_fft.h"
#include <math.h>
int fft(complex* x)
{
int i =0;
int ip = 0;
int j = 0;
int jm1 = 0;
int k = 0;
int l = 0;
int le = 0;
int le2 = 0;
int nm1 = FFT_N-1;
int nd2 = FFT_N/2;
int m = (int)(log(FFT_N)/log(2));
complex t = {.real = 0.0, .imag = 0.0};
complex u = {.real = 1.0, .imag = 0.0};
complex s = {.real =0.0, .imag = 0.0};
/*bit reversal sorting of x_in*/
j = nd2;
for(i = 1; i < nm1; i ++)
{
if(i < j) // while the i < j , transverse x,x,
{
t=x;
x = x;
x = t;
}
//Rather mechod to bit reverse arrangement
k = nd2;
while(k <= j)
{
j = j - k; //set the big end bit to 0
k = k/2; //the next bit , from up bit k to next bit k/2...
}
j = j + k; //set the big end bit to 1, as k > j
}
/******************************************************
*
* DFT routine
*
*******************************************************/
for(l =1 ; l <= m; l ++)
{
le = (int)(1<<l);
le2 = le/2;
u.real = 1.0;
u.imag = 0.0;
s.real = cos(PI/le2);
s.imag = -sin(PI/le2);
for(j = 1 ; j <= le2; j++)//for each sub DFT
{
jm1 = j -1;
for(i = jm1; i <= nm1; i = i + le) //loop for each Butterfly
{
ip = i + le2;
t.real = x.real * u.real - x.imag * u.imag;
t.imag = x.real * u.imag + x.imag * u.real;
x.real = x.real - t.real;
x.imag = x.imag - t.imag;
x.real = x.real + t.real;
x.imag = x.imag + t.imag;
}
t.real = u.real;
u.real = t.real * s.real - u.imag * s.imag;
u.imag = t.real * s.imag + u.imag * s.real;
}
}
return 0;
}
/********************************************************************************
*
* INVERSE FAST FOURIER TRANSFORM FUNCTION
*
********************************************************************************/
int ifft(complex * x)
{
int i = 0;
int k = 0;
int n = FFT_N ;
int nm1 = FFT_N - 1;
//change the sign of imag part of x[]
for( k = 0; k < n; k++)
{
x.imag = -x.imag;
}
//forward FFT
fft(x);
for(i = 0 ; i< n; i++) //divide the time domain by N andchange the sign of x[].imag.
{
x.real = x.real/n;
x.imag = -x.imag/n;
}
return 0;
}
4.测试样例：
5. 测试的效果：
6.致谢： 感谢南京沁恒公司对我的厚爱，和沁恒的缘分很长，2014年做Arduino的板子，使用的是CH340G，之前也申请过CH32F103的开发板，沁恒也是支持。这次的比赛我也很幸运得到了开发板的物质支持。向沁恒的工程师们致敬！谢谢！
这个运算速度怎么？要几个mS GSDDDD 发表于 2021-8-30 21:56
这个运算速度怎么？要几个mS
测试的是4096个点的。这个库严格讲仅仅是实现了FFT算法，并没有进行优化，偏向原理的学习。但是对于一般的使用应该是可以的。运行的时间我没有测试。使用计时器计时，我还需要学习。 catnull 发表于 2021-8-30 23:47
测试的是4096个点的。这个库严格讲仅仅是实现了FFT算法，并没有进行优化，偏向原理的学习。但是对于一般 ...
运行时间，用一个IO指示即可。FFT，MCU都可以计算，关键是速度。 计算的速度是多少？ 这个是官网的库函数吗？ 支持一下。 这个是最近才申请的吗 fft的计算很浪费资源。 快速傅里叶和离散傅里叶一样吗 72Mhz的计算量真是一般呢 做个数据解析吧。 能把matlab的代码转换过来吗 QT有显示吗 这个芯片有硬件计算吗 内部的ram够不够计算呢 代码实现是一致的吗 ch32v103最快好像是80Mhz的频率。 ch32v307可以到180Mhz的主频。 怎么使用lib的函数呢