前言

最近工作移植PPG算法，将MATLAB上代码移植到嵌入式设备上去。因为心率算法利用FFT实现会较为简单，所以又重新了解了一下大学里学的FFT，并写了C语言实现MATLAB的FFT接口的代码。看了好多都是用的递归写的，这样对于算法复杂度来说还是挺大的，这里参考了这篇大佬的文章，将大佬的代码稍加修改，整体效果还是不错的。

FFT介绍

1.DFT与FFT

DFT一般是指离散傅里叶变换（Discrete Fourier Transform，DFT）是信号分析的最基本方法，傅里叶变换是傅里叶分析的核心，通过它把信号从时间域变换到频率域，进而研究信号的频谱结构和变化规律。

离散傅里叶变换（DFT），是傅里叶变换在时域和频域上都呈现离散的形式，将时域信号的采样变换为在离散时间傅里叶变换（DTFT）频域的采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作经过周期延拓成为周期信号再作变换。在实际应用中通常采用快速傅里叶变换以高效计算DFT。 ——百度

FFT一般指快速傅里叶变换 (fast Fourier transform), 即利用计算机计算离散傅里叶变换（DFT)的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。——百度

这里可以这么理解，FFT就是计算机加速DFT计算过程的算法的称呼。

2.FFT算法介绍

DFT公式

\[X(k) = DFT[X(n)] = \sum_{n = 0}^{N-1} x(n)W_N^{nk} \]

旋转因子

计算DFT时，如采用暴力计算，则：

\[f(x)=a_0+a_1x+a_2x^2+a_3x^3+...+a_nx^n \]

这里我们可以引入一些特别的x，让他们的若干次平方之后结果为1。

首先我们可以想到1和-1满足条件，1的多少次方都是1，-1的偶数次方是1，但是又出现了一个新的问题，那就是1和-1只有两个数，而多项式需要的是n个不同的数。

然后我们想到了复数，1，i，-1，-i，但是这也只有4个数，依然无法满足条件。

最后傅里叶引入了一个单位圆。

然后把园N等分后（N是2的整数次幂）就得到了

\[W_N^k=cos(-\frac{2\pi k}{N})+isin(\frac{2\pi k}{N}) \]

其中k就是圆的第k等份

由旋转因子的对称性、均匀性、周期性，我们可以整理出旋转因子的两个很重要的性质。

\[W_N^k=W_{2N} ^{2k} \]

\[W_N^{k+\frac{\pi}{2}}=-W_{N} ^{k} \]

将旋转因子带入计算

具体推理过程可以参考大佬的文章

得到两个重要公式：

\[X(W_N^k)=F(W_\frac N 2 ^k)+W_N^k G(W_\frac N 2 ^k) \]

\[X(W_N^{k+\frac{\pi}{2}})=F(W_\frac N 2 ^k)-W_N^k G(W_\frac N 2 ^k) \]

蝶形运算

先以一个4点蝶形运算图为例说明：详细可以参考这篇文章

第一次了解到这个图的时候是在大三学数字信号与系统，当时期末考试每年必有画这个图的题目，当时也不懂这个图只知道如何画，刷了好多道题，题目到会做了，但这个图所带来的意义一点也不明白。最近深入了解了后才明白这个图的意义，也是代码的核心部分。

当k=0时

\[\begin{aligned} X(W_4^0)&=F(W_2^0)+W_4^0 G(W_2^0)\\ &=\sum_{n = 0}^{1}x(2n)(W_2^0)^n+W_4^0\sum_{n = 0}^{1}x(2n+1)(W_2^0)^n \end{aligned} \]

\[\begin{aligned} F(W_2^0) &=F_f(W_1^0)+W_2^0G_f(W_1^0) \\ &=\sum_{n = 0}^{0}x(4n)(W_1^0)^n+W_2^0\sum_{n = 0}^{1}x(4n+2)(W_1^0)^n\\ &=x(0)+x(2)W_2^0 \end{aligned} \]

\[\begin{aligned} G(W_2^0) &=F_g(W_1^0)+W_2^0G_g(W_1^0) \\ &=\sum_{n = 0}^{0}x(4n+1)(W_1^0)^n+W_2^0\sum_{n = 0}^{1}x(4n+3)(W_1^0)^n\\ &=x(1)+x(3)W_2^0 \end{aligned} \]

其中 \(F_f\) , \(G_f\)为\(F\)的递归函数， \(F_g\) , \(G_g\)为\(G\)的递归函数。

所以可得\(X(W^0_ {4})\)为:

\[x(W_4^0)=x(0)+X(2)W_2^0+W_4^0(x(1)+x(3)W_2^0) \]

这里就为FFT递归版本算法。网上有很多递归版本的写法，也比较容易，这里采用一种迭代\(for\)循环的方式。

注意观察前后位置的变化\(x[0]\quad x[2] \quad x[1] \quad x[3]\)变成\(X[0]\quad X[1] \quad X[2] \quad X[3]\)再看一手8点的。

注意观察：

它们分割后排序的二进制表示刚好就是没分割前的二进制表示的相反结果

到此，算法基本讲解完成，接下来展示代码。

C语言代码

1.定义复数结构体

 /* 复数结构体 */
typedef struct 
{
    double real;
	double imag;
}Complex;

2.定义复数运算

/* 复数加法 */
void Complex_ADD(Complex a, Complex b, Complex *c)
{
	c->real = a.real + b.real; 
	c->imag = a.imag + b.imag;
}
/* 复数减法 */
void Complex_SUB(Complex a, Complex b, Complex *c)
{
	c->real = a.real - b.real;
	c->imag = a.imag - b.imag;
}
/* 复数乘法 */
void Complex_MUL(Complex a, Complex b, Complex *c)
{
	c->real = a.real * b.real - a.imag * b.imag;
	c->imag = a.real * b.imag + a.imag * b.real;
}
/* 复数绝对值 */
double Complex_ABS(Complex *a)
{
	double b;
	b = sqrt((a->real)*(a->real)+(a->imag)*(a->imag));
	return b;
}

3.FFT主函数

static double PI=4.0*atan(1); //定义π 因为tan(π/4)=1 所以arctan（1）*4=π，增加π的精度