快速傅里叶变换 FFT C++

#include <math.h>    
#define DOUBLE_PI   6.283185307179586476925286766559    
   
// 快速傅里叶变换    
// data 长度为 (2 * 2^n), data 的偶位为实数部分, data 的奇位为虚数部分    
// isInverse表示是否为逆变换    
void FFT(double * data, int n, bool isInverse = false)   
{   
    int mmax, m, j, step, i;   
    double temp;   
    double theta, sin_htheta, sin_theta, pwr, wr, wi, tempr, tempi;   
    n = 2 * (1 < n);   
    int nn = n >> 1;   
    // 长度为1的傅里叶变换, 位置交换过程    
    j = 1;   
    for(i = 1; i  n; i += 2)   
    {   
        if(j > i)   
        {   
            temp = data[j - 1];   
            data[j - 1] = data[i - 1];   
            data[i - 1] = temp;   
            data[j] = temp;   
            data[j] = data[i];   
            data[i] = temp;   
        }   
        // 相反的二进制加法    
        m = nn;   
        while(m >= 2 && j > m)   
        {   
            j -= m;   
            m >>= 1;   
        }   
        j += m;   
    }   
    // Danielson - Lanczos 引理应用    
    mmax = 2;   
    while(n > mmax)   
    {   
        step = mmax < 1;   
        theta = DOUBLE_PI / mmax;   
        if(isInverse)   
        {   
            theta = -theta;   
        }   
        sin_htheta = sin(0.5 * theta);   
        sin_theta = sin(theta);   
        pwr = -2.0 * sin_htheta * sin_htheta;   
        wr = 1.0;   
        wi = 0.0;   
        for(m = 1; m  mmax; m += 2)   
        {   
            for(i = m; i = n; i += step)   
            {   
                j = i + mmax;   
                tempr = wr * data[j - 1] - wi * data[j];   
                tempi = wr * data[j] + wi * data[j - 1];   
                data[j - 1] = data[i - 1] - tempr;   
                data[j] = data[i] - tempi;   
                data[i - 1] += tempr;   
                data[i] += tempi;   
            }   
            sin_htheta = wr;   
            wr = sin_htheta * pwr - wi * sin_theta + wr;   
            wi = wi * pwr + sin_htheta * sin_theta + wi;   
        }   
        mmax = step;   
    }   
}   
   
  /*  输入数据为data，data是一组复数，偶数位存储的是复数的实数部分，奇数位存储的是复数的虚数部分。data的长度与n相匹配。注意：这里的n并非是data的长度，data的实际长度为(2 * 2^n),存储了N = 2^n个复数。  
  
    输出也存放在data中。  
  
    以正向傅里叶变换为例，作为输入data中存储的是以delta为时间间隔时域函数的振幅抽样值。经过函数计算后data中存放输出，存储的是以1/(N * delta)为频率间隔频域像函数值。频率范围为0Hz,1/(N * delta),2/(N * delta) ... (N / 2 - 1) / N * delta, +/- 1 / delta, -(N / 2 - 1) / N * delta ... -2/(N * delta), -1/(N * delta)。注意这是一个中间大两边小的排列。  
  
    如果将isInverse设置为true则计算逆傅里叶变换。*/