前言

最近工作移植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=π，增加π的精度

/**
 * @description: FFT演算法（非遞迴）
 * @param {cs_uint32*} signal_in: 需要fft訊號輸入
 * @param {cs_uint32} signal_len: 輸入訊號長度
 * @param {Complex*} fft_out: fft完成後輸出結果
 * @param {cs_uint32} fft_point:  取樣點數
 * @return {*} 
 */
uint32_t FFT(uint32_t* signal_in, uint32_t signal_len, Complex* fft_out, uint32_t fft_point)
{
	Complex *W;//旋轉因子Wn
	uint32_t i,j,k,m;
	Complex temp1,temp2,temp3;//用於交換中間量
	double series;//序列級數
	if(signal_len < fft_point)//取樣點數與訊號資料對比 
	{
		for(i=0;i < signal_len;i++)
		{
			fft_out[i].real = signal_in[i];
			fft_out[i].imag = 0;
		}
		for(i=signal_len;i<fft_point;i++)//補0
		{
			fft_out[i].real = 0;
			fft_out[i].imag = 0;
		}
	}else if(signal_len == fft_point)
	{
		for(i=0;i < signal_len;i++)
		{
			fft_out[i].real = signal_in[i];
			fft_out[i].imag = 0;
		}
	}
	else
	{
		for(i=0;i < fft_point;i++)
		{
			fft_out[i].real = signal_in[i];
			fft_out[i].imag = 0;
		}
	}
	if(fft_point%2 != 0)
	{
		return 1;
	}
	
	W = (Complex *)malloc(sizeof(Complex) * fft_point);
	if (W == NULL)
	{
		return 1;
	}
	
	for (i = 0; i < fft_point; i++)
	{
		W[i].real = cos(2*PI/fft_point*i);	 //尤拉表示的實部
		W[i].imag = -1*sin(2*PI/fft_point*i);//尤拉表示的虛部
	}
	
	for(i = 0;i < fft_point;i++)
	{
		k = i;
		j = 0;
		series = log(fft_point)/log(2); //算出序列的級數
		while((series--) > 0)//利用按位與以及迴圈實現碼位顛倒
		{
			j = j << 1;
			j|=(k & 1);
			k = k >> 1;
		}
		if(j > i) //將x(n)的碼位互換
		{
			temp1 = fft_out[i];
			fft_out[i] = fft_out[j];
			fft_out[j] = temp1;
		}
	}
	series = log(fft_point)/log(2); //蝶形運算的級數
	
	for(i = 0; i<series;i++)
	{
		m = 1<<i; //移位 每次都是2的指數的形式增加，其實也可以用m=2^i代替
		for(j = 0;j<fft_point;j+=2*m) //一組蝶形運算，每一組的蝶形因子乘數不同
		{
			for(k=0;k<m;k++)//蝶形結點的距離  一個蝶形運算 每個組內的蝶形運算
			{
				Complex_MUL(fft_out[k+j+m],W[fft_point*k/2/m],&temp1);
				Complex_ADD(fft_out[j+k],temp1,&temp2);
				Complex_SUB(fft_out[j+k],temp1,&temp3);
				fft_out[j+k] = temp2;
				fft_out[j+k+m] = temp3;
			}
		}
	}
	free(W);
	return 0;
}