fft.h

基于快速傅立叶变换的c语言程序

#ifndef _FFT_HERDER
#define _FFT_HERDER

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <complex>
using std::complex;
using std::polar;

// 常用宏
#define SQ2(a)      ((a)*(a))
#define ABS(a)      ((a)<0 ?-(a):(b))
#define MIN(a, b)   ((a)<(b)?(a):(b))
#define MAX(a, b)   ((a)>(b)?(a):(b))

typedef complex<double> complexD;       // 实虚部是double型的复数

static const double PI=3.1415926;
static const int N = 1024;              // 计算点阵
static int *rev;                        // 码位倒置数组
static complexD *w_fact, *w_fact_r;     // 正反FFT的W旋转因子

/*---------------------------------------------------------------------------------------------
 *                              FFT算法初始化函数
 *  功能：
 *      加速起见，FFT需要的3个辅助数组：rev，w_fact，w_fact_r。因此请务必在main函数的
 *      起始就调用FFT_INIT。程序终结处，调用FFT_FREE，来释放这3个动态数组。
 *  说明：
 *      对于每次不同的n长度的傅里叶变换，那么就要FFT_INIT和FFT_FREE夹住之。因为FFT_INT
 *      计算码位倒置数组和W旋转因子，是针对某个n值的。
 *--------------------------------------------------------------------------------------------*/
void FFT_INIT(int n)
{
    int i, j, k, p, q;
    for (k=0, j=1; j < n; j<<=1, k++);
    if (j != n) {
        fprintf(stderr, "%d must be power of 2!\n", n);
        exit(1);
    }
    
    rev = new int [n];
    w_fact = new complexD [n];
    w_fact_r = new complexD [n];

    rev[0] = 0;
    for (i = 1, p = 0; i < n; i++)
    {
        for (q=1<<(k-1); q; q>>=1)
        {
            if (q & p) {
                p &= ~q;
            } else {
                p |= q;
                break;
            }
        }
        rev[i] = p;
    }

    j = n >> 1;
    for (i = 0; i < j; i++)
    {
        w_fact[i] = polar(1.0, -2*i*PI/n);
        w_fact_r[i] = polar(1.0, 2*i*PI/n);
    }
}
void FFT_FREE()
{
    delete [] rev;
    delete [] w_fact;
    delete [] w_fact_r;
}

/*---------------------------------------------------------------------------------------------
 *                          1维和2维相关的FFT函数
 *  说明：
 *      _FFT:   纯数学的1维FFT，相当于Matlab的fft
 *      FFT:    带中心基频移动的FFT，正向变换相当于Matlab的：fft + fftshift；
 *              逆向变换，相当于 fftshift + fft
 *      _FFT2D: 纯数学的2维FFT，相当于Matlab的fft2
 *      FFT2D:  带中心基频移动的FFT，正向变换相当于Matlab的：fft2 + fftshift；
 *              逆向变换，相当于 fftshift + fft2
 *      FFT2D:  不同参数表，相当于重载，表示变换的输出矩阵不是原输入矩阵。
 *  flag参数：1表示正向，-1表示逆向
 *--------------------------------------------------------------------------------------------*/
void _FFT(complexD *A, int n, int flag) 
{
    int i, j, k, l, p, q;
    int M, Nv2 = n/2;
    int step1, step2, stp;
    complexD temp, W, *w_ptr;
    flag = (flag<0 ? -1 : 1);
    w_ptr= (flag<0 ? w_fact_r : w_fact);
    for (M = 0; 1<<M < n; ++M);

    for (i = 0; i < n; i++)
    {
        j = rev[i];
        if (i < j) {
            temp = A[i];
            A[i] = A[j];
            A[j] = temp;
        }
    }
    for (i = 1; i <= M; i++)
    {
        step2 = 1 << i;
        step1 = step2 >> 1;
        stp = 1 << (M-i);
        l = 0;
        W = w_ptr[0];

        for (j = 0; j < step1; j++, W=w_ptr[l+=stp])
            for (k = 0; k < n; k += step2)
            {
                p = k + j;
                q = p + step1;
                temp = W * A[q];
                A[q] = A[p] - temp;
                A[p] += temp;
            }
    }

    double ceof = 1./n;
    if (flag < 0)
        for(i = 0; i < n; ++i) A[i] *= ceof;
}
void FFT(complexD *A, int n, int flag)  // 带有中心基频移动的FFT 
{
    int i, j;
    int ct = n >> 1;
    complexD t;

    if (flag < 0)
        for(i=0, j=ct; i < ct; ++i, ++j)
        { t = A[i]; A[i] = A[j]; A[j] = t; }

    _FFT(A, n, flag);

    if (flag > 0)
        for (i=0, j=ct; i < ct; ++i, ++j)
        { t = A[i]; A[i] = A[j]; A[j] = t; }
}

void _FFT2D(complexD *A[N], int flag)
{
    int i, j;
    complexD t;

    for(i = 0; i < N; i++)
        _FFT(A[i], N, flag);
    for(i = 1; i < N; i++)  // In-Place 矩阵转置
        for(j = 0; j < i; j++)
        { t = A[i][j]; A[i][j] = A[j][i]; A[j][i] = t; }
    for(i = 0; i < N; i++)
        _FFT(A[i], N, flag);
    for(i = 1; i < N; i++)
        for(j = 0; j < i; j++)
        { t = A[i][j]; A[i][j] = A[j][i]; A[j][i] = t; }
}
void FFT2D(complexD *A[N], int flag)
{
    int i, j;
    complexD t;

    for(i = 0; i < N; i++)
        FFT(A[i], N, flag);
    for(i = 1; i < N; i++)  // In-Place 矩阵转置
        for(j = 0; j < i; j++)
        { t = A[i][j]; A[i][j] = A[j][i]; A[j][i] = t; }
    for(i = 0; i < N; i++)
        FFT(A[i], N, flag);
    for(i = 1; i < N; i++)
        for(j = 0; j < i; j++)
        { t = A[i][j]; A[i][j] = A[j][i]; A[j][i] = t; }
}
void FFT2D(complexD *B[N], complexD *A[N], int flag)  // 重载 B = FFT2D(A)
{
    int i, j;
    complexD t;

    for(i = 0; i < N; i++)
        for(j = 0; j < N; j++)
            B[i][j] = A[j][i];
    for(i = 0; i < N; i++)
        FFT(B[i], N, flag);
    for(i = 1; i < N; i++)
        for(j = 0; j < i; j++)
        { t = B[i][j]; B[i][j] = B[j][i]; B[j][i] = t; }
    for(i = 0; i < N; i++)
        FFT(B[i], N, flag);
}

#endif