快速傅里叶变换(FFT)算法C++实现代码

快速傅里叶变换(FFT)算法C++实现代码

#include

#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;

data = 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 - tempi;

data[i - 1] += tempr;

data += 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 / d elta, -(N / 2 - 1) / N * delta ... -2/(N * delta), -1/(N * delta)。注意这是一个中间大两边小的排列。

如果将isInverse设置为true则计算逆傅里叶变换。