FFT算法的完整DSP實現

傅里葉變換或者FFT的理論參考：

[1] http://www.dspguide.com/ch12/2.htm

      The Scientist and Engineer's Guide to Digital Signal Processing,   By Steven W. Smith, Ph.D.

[2] http://blog.csdn.net/v_JULY_v/article/details/6196862，可當作[1]的中文參考

[3] 任意一本數字信號處理教材，上面都有詳細的推導DCT求解轉換為FFT求解的過程

[4] TI文檔：基於TMS320C64x+DSP的FFT實現。 使用baidu/google可以搜索到。

1. 有關FFT理論的一點小小解釋

關於FFT這里只想提到兩點：

（1）DFT變換對的表達式（必須記住）

          —— 稱旋轉因子

（2）FFT用途——目標只有一個，加速DFT的計算效率。

DFT計算X(k)需要N^2次復數乘法和N(N-1)次復數加法；FFT將N^2的計算量降為。

“FFT其實是很難的東西，即使常年在這個領域下打拼的科學家也未必能很好的寫出FFT的算法。”——摘自參考上面提供的參考文獻[1]

因此，我們不必太過糾結於細節，當明白FFT理論后，將已有的算法挪過來用就OK了，不必為閉着教材寫不出FFT而郁悶不堪。

FFT的BASIC程序偽代碼如下：

IFFT的BASIC程序偽代碼如下（IFFT通過調用FFT計算）：

FFT算法的流程圖如下圖，總結為3過程3循環：

（1）3過程：單點時域分解（倒位序過程） + 單點時域計算單點頻譜 + 頻域合成

（2）3循環：外循環——分解次數，中循環——sub-DFT運算，內循環——2點蝶形算法

分解過程或者說倒位序的獲得參考下圖理解：

 

2. FFT的DSP實現

 

下面為本人使用C語言實現的FFT及IFFT算法實例，能計算任意以2為對數底的采樣點數的FFT，算法參考上面給的流程圖。

 

/*
 * zx_fft.h
 *
 *  Created on: 2013-8-5
 *      Author: monkeyzx
 */

#ifndef ZX_FFT_H_
#define ZX_FFT_H_

typedef float          FFT_TYPE;

#ifndef PI
#define PI             (3.14159265f)
#endif

typedef struct complex_st {
	FFT_TYPE real;
	FFT_TYPE img;
} complex;

int fft(complex *x, int N);
int ifft(complex *x, int N);
void zx_fft(void);

#endif /* ZX_FFT_H_ */

 

/*
 * zx_fft.c
 *
 * Implementation of Fast Fourier Transform(FFT)
 * and reversal Fast Fourier Transform(IFFT)
 *
 *  Created on: 2013-8-5
 *      Author: monkeyzx
 */

#include "zx_fft.h"
#include <math.h>
#include <stdlib.h>

/*
 * Bit Reverse
 * === Input ===
 * x : complex numbers
 * n : nodes of FFT. @N should be power of 2, that is 2^(*)
 * l : count by bit of binary format, @l=CEIL{log2(n)}
 * === Output ===
 * r : results after reversed.
 * Note: I use a local variable @temp that result @r can be set
 * to @x and won't overlap.
 */
static void BitReverse(complex *x, complex *r, int n, int l)
{
	int i = 0;
	int j = 0;
	short stk = 0;
	static complex *temp = 0;

	temp = (complex *)malloc(sizeof(complex) * n);
	if (!temp) {
		return;
	}

	for(i=0; i<n; i++) {
		stk = 0;
		j = 0;
		do {
			stk |= (i>>(j++)) & 0x01;
			if(j<l)
			{
				stk <<= 1;
			}
		}while(j<l);

		if(stk < n) {             /* 滿足倒位序輸出 */
			temp[stk] = x[i];
		}
	}
	/* copy @temp to @r */
	for (i=0; i<n; i++) {
		r[i] = temp[i];
	}
	free(temp);
}

/*
 * FFT Algorithm
 * === Inputs ===
 * x : complex numbers
 * N : nodes of FFT. @N should be power of 2, that is 2^(*)
 * === Output ===
 * the @x contains the result of FFT algorithm, so the original data
 * in @x is destroyed, please store them before using FFT.
 */
int fft(complex *x, int N)
{
	int i,j,l,ip;
	static int M = 0;
	static int le,le2;
	static FFT_TYPE sR,sI,tR,tI,uR,uI;

	M = (int)(log(N) / log(2));

	/*
	 * bit reversal sorting
	 */
	BitReverse(x,x,N,M);

	/*
	 * For Loops
	 */
	for (l=1; l<=M; l++) {   /* loop for ceil{log2(N)} */
		le = (int)pow(2,l);
		le2 = (int)(le / 2);
		uR = 1;
		uI = 0;
		sR = cos(PI / le2);
		sI = -sin(PI / le2);
		for (j=1; j<=le2; j++) {   /* loop for each sub DFT */
			//jm1 = j - 1;
			for (i=j-1; i<=N-1; i+=le) {  /* loop for each butterfly */
				ip = i + le2;
				tR = x[ip].real * uR - x[ip].img * uI;
				tI = x[ip].real * uI + x[ip].img * uR;
				x[ip].real = x[i].real - tR;
				x[ip].img = x[i].img - tI;
				x[i].real += tR;
				x[i].img += tI;
			}  /* Next i */
			tR = uR;
			uR = tR * sR - uI * sI;
			uI = tR * sI + uI *sR;
		} /* Next j */
	} /* Next l */

	return 0;
}

/*
 * Inverse FFT Algorithm
 * === Inputs ===
 * x : complex numbers
 * N : nodes of FFT. @N should be power of 2, that is 2^(*)
 * === Output ===
 * the @x contains the result of FFT algorithm, so the original data
 * in @x is destroyed, please store them before using FFT.
 */
int ifft(complex *x, int N)
{
	int k = 0;

	for (k=0; k<=N-1; k++) {
		x[k].img = -x[k].img;
	}

	fft(x, N);    /* using FFT */

	for (k=0; k<=N-1; k++) {
		x[k].real = x[k].real / N;
		x[k].img = -x[k].img / N;
	}

	return 0;
}

/*
 * Code below is an example of using FFT and IFFT.
 */
#define  SAMPLE_NODES              (128)
complex x[SAMPLE_NODES];
int INPUT[SAMPLE_NODES];
int OUTPUT[SAMPLE_NODES];

static void MakeInput()
{
	int i;

	for ( i=0;i<SAMPLE_NODES;i++ )
	{
		x[i].real = sin(PI*2*i/SAMPLE_NODES);
		x[i].img = 0.0f;
		INPUT[i]=sin(PI*2*i/SAMPLE_NODES)*1024;
	}
}

static void MakeOutput()
{
	int i;

	for ( i=0;i<SAMPLE_NODES;i++ )
	{
		OUTPUT[i] = sqrt(x[i].real*x[i].real + x[i].img*x[i].img)*1024;
	}
}

void zx_fft(void)
{
	MakeInput();

	fft(x,128);
	MakeOutput();

	ifft(x,128);
	MakeOutput();
}

程序在TMS320C6713上實驗，主函數中調用zx_fft()函數即可。

 

FFT的采樣點數為128，輸入信號的實數域為正弦信號，虛數域為0，數據精度定義FFT_TYPE為float類型，MakeInput和MakeOutput函數分別用於產生輸入數據INPUT和輸出數據OUTPUT的函數，便於使用CCS 的Graph功能繪制波形圖。這里調試時使用CCS v5中的Tools -> Graph功能得到下面的波形圖（怎么用自己琢磨，不會的使用CCS 的Help）。

輸入波形

輸入信號的頻域幅值表示

FFT運算結果

對FFT運算結果逆變換(IFFT)

如何檢驗運算結果是否正確呢？有幾種方法：

（1）使用matlab驗證，下面為相同情況的matlab圖形驗證代碼

 

SAMPLE_NODES = 128;
i = 1:SAMPLE_NODES;
x = sin(pi*2*i / SAMPLE_NODES);

subplot(2,2,1); plot(x);title('Inputs');
axis([0 128 -1 1]);

y = fft(x, SAMPLE_NODES);
subplot(2,2,2); plot(abs(y));title('FFT');
axis([0 128 0 80]);

z = ifft(y, SAMPLE_NODES);
subplot(2,2,3); plot(abs(z));title('IFFT');
axis([0 128 0 1]);

 

（2）使用IFFT驗證：輸入信號的FFT獲得的信號再IFFT，則的到的信號與原信號相同

可能大家發現輸入信號上面的最后IFFT的信號似乎不同，這是因為FFT和IFFT存在精度截斷誤差（也叫數據截斷噪聲，意思就是說，我們使用的float數據類型數據位數有限，沒法完全保留原始信號的信息）。因此，IFFT之后是復數（數據截斷噪聲引入了虛數域，只不過值很小），所以在繪圖時使用了計算幅值的方法，

C代碼中：

OUTPUT[i] = sqrt(x[i].real*x[i].real + x[i].img*x[i].img)*1024;

 

matlab代碼中：

 

subplot(2,2,3); plot(abs(z));title('IFFT');

 

所以IFFT的結果將sin函數的負y軸數據翻到了正y軸。另外，在CCS v5的圖形中我們將顯示信號的幅度放大了1024倍便於觀察，而matlab中沒有放大。