C语言实现FFT算法 - 开发手记

C语言实现FFT算法 - 开发手记

Realize FFT transformation by C language.

by Shun Fu,
DSP Lab of Northeastern University.

In order to do frequency analyze for an input signal, we need to do the DFT transform for the data sampled from our codec(AD/DA).



Read more about DFT: http://en.wikipedia.org/wiki/Discrete_Fourier_transform
FFT is a more efficiency algorithm for calculating the DFT.

Here I don't describe what is FFT algorithm, you can check it on wikipedia.com

Here is a figure explaining the principle of 3 order （N = 2^3 = 8points） FFT algorithm.





Figure 1: Picked as time, 8 point, FFT signal stream figure.

Actually, for the first column of x[i], they are position exchanged by the column: x[0] , x[1],...,x[7].

So, we need to firstly change the order of our sample data.

Step 1: Order your date sampled from you ADC:
In case of 3 order, 8 datas: The mechanism of ordering is:
Assume x[i] is the array of our original data from ADC, and x_0[i] is the array that was ordered by us.
x_0[0] = x_0[000] <= x[000] = x[0]
x_0[1] = x_0[001] <= x[100] = x[4]
x_0[2] = x_0[010] <= x[010] = x[2]
x_0[3] = x_0[011] <= x[110] = x[6]
x_0[4] = x_0[100] <= x[001] = x[1]
x_0[5] = x_0[101] <= x[101] = x[5]
x_0[6] = x_0[110] <= x[011] = x[3]
x_0[7] = x_0[111] <= x[111] = x[7]
if order equals 4, x_0[12] = x_0[1100] <= x[0011] = x[3].
reference:
Order data using the algorithm of Interchange Function of Computer Architecture
用“计算机结构”的“互联函数”思想构造FFT算法中每层之间的数据地址构造问题

C Code to realize this order exchange mechanism:



In above case, we can call the function by: exch_bit(12,4) , and we can get 3 which is the return value.

The input parameter 12 is the number of original data and 4 is the order. 3 is the new order number that we need. In case of order = 4, we can use the instruction below to exchange the order of original data:
order = 4;
N = xp(2,4);
for(j = 0, j < N; j++ )
{
x_0[j] = x[exch_bit()]；
}

Step 2：Code for FFT calculation:

The factor WN^kn contains Re part and Im part, the code for calculate the FFT split it into Re and Im. It calculate them separately. r below is k*n.

Please see the C code for pool of FFT calculation below, you may need to refer to the figure1 above.



Explain for the code above：
i. r is for calculate factor WN^kn, r = kn. In figure 1, each array of x_i[j] is the result of x_i-1[xx] multiple the factor WN^r and add x_i-1[xxx].

3 order, i = 1, j = 0, => r=0
3 order, i = 1, j = 1, => r=4
3 order, i = 2, j = 0, => r=0
3 order, i = 2, j = 1, => r=2
3 order, i = 2, j = 2, => r=4
3 order, i = 2, j = 3, => r=6
ii. cpx_x_Re is a function that calculate the multiple operation of two complex number, and return the Re part of the result. for example: For two complex number: c_1 = (a,b), c_2 =
(c,d).c_3 = c_1 * c_2 = (ac - bd) + i(ad + bc) = [(ac - bd), (ad + bc)]
Calling of cpx_x_Re and cpx_x_Im can do the operation above:

cpx_x_Re(a,b,c,d) = ac - bd
cpx_x_Im(a,b,c,d) = ad + bc


FFT transformation can be realized by different ways. Here is the way chosing 2 as the base as described by the classical signals and systems book written by Alan V. Oppenheim, and pick out data as time. Most of
the MP3 WMA codec uses the FFT transformation of picking out as frequency which has more efficiency. And some of the FFT transform do calculate choosing 4 as the base and do base-2 again. they can decrease N/2 multiple operations actually.

Reference:
http://blog.csdn.net/xcgspring/article/details/4749075