SOFTIST 编程方法筆記 →目録

快速傅氏变换(算法)

Fast Fourier Transform, FFT

快速傅立叶变换。

离散傅立叶变换 (DFT)的定义如下。由于级数的计算量过大，花费大量CPU处理时间，在实时处理领域里并没有多少实用价值。

    　  　N-1　　   jk                -2πi/N    
    A  =  Σ   a   W   ，      W  =  e   　         
    　k 　j=0　 j   N  　   　  N 　 　   　         

    　  　 1  N-1        -jk    
    a  =  --- Σ   A   W       
     k 　  N  j=0   j   N      

一般Cooley-Tukey型FFT算法实现的快速傅立叶变换被广泛使用。不过使用这个算法有一个限制，就是参加运算的离散数字的个数必须为２的整数倍。下面是Cooley-Tukey型FFT的C程序的関数。当然逆変换（IFFT）时也可以使用这个函数。n是数据的个数，２的整数倍。ar[], ai[]是复数的实部和虚部。 theta是运算系数，其值为 -8 * atan(1.0) / n，逆変換時它的值应设成8  * atan(1.0) / n。

void fft(int n, double theta, double ar[], double ai[])
{
    int m, mh, i, j, k;
    double wr, wi, xr, xi;
    for (m = n; (mh = m >> 1) >= 1; m = mh) 
    {
        for (i = 0; i < mh; i++) 
        {
            wr = cos(theta * i);
            wi = sin(theta * i);
            for (j = i; j < n; j += m) 
            {
                k = j + mh;
                xr = ar[j] - ar[k];
                xi = ai[j] - ai[k];
                ar[j] += ar[k];
                ai[j] += ai[k];
                ar[k] = wr * xr - wi * xi;
                ai[k] = wr * xi + wi * xr;
            }
        }
        theta *= 2;
    }
    
    i = 0;
    for (j = 1; j < n - 1; j++) 
    {
        for (k = n >> 1; k > (i ^= k); k >>= 1);
        if (j < i) 
        {
            xr = ar[j];
            xi = ai[j];
            ar[j] = ar[i];
            ai[j] = ai[i];
            ar[i] = xr;
            ai[i] = xi;
        }
    }
}

使用例：低通濾波器 (Low-pass filter, LPF)

#define FFT_SIZE 1024
void LPF(ar, ai)    
{
    double theta = -8 * atan(1.0) / FFT_SIZE;
    fft(FFT_SIZE, theta, ar, ai);    

    //用零值替代高频部分，以实现低通滤波
    for (i = iFreqHightIndex; i < FFT_SIZE / 2; i ++)
    {
        ar[i] = ai[i] = ar[FFT_SIZE - i - 1] = ai[FFT_SIZE - i - 1] = 0.0;
    }

    theta = 8 * atan(1.0) / FFT_SIZE;
    fft(FFT_SIZE, theta, ar, ai);
}