这次接着上次说, 花了点时间把按时间抽取基-2的FFT思路用C++代码实现了下. 可以利用递归的思路来实现.

而且, 就像之前理解的FFT那样, 直接在一个序列的存储空间中进行操作, 这里存在vector中. 要实现这个特性, 调用时必须传入当前处理的那一段序列的序号范围, 所以第一次调用FFT时又套了个壳(myFFT). 有遗憾的是没有实现旋转因子的复用, 也就是重复出现的旋转因子都进行了重复计算. 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
#include "mylib\mycomplex.h"
#include <iostream>
#include <vector>
#include <cmath>

using namespace std;
using namespace mylib;

vector<Complex> FFT(vector<Complex> &xn, int begin, int end)
{
	int length = end - begin + 1;//计算该段序列长度
	if (length > 2)//说明需要分组
	{
		for (int i = begin + 1; i <= begin + length / 2; i++)
		{
			xn.insert(xn.end() - (xn.size() - end - 1), xn[i]);
			xn.erase(xn.begin() + i);
		}
		FFT(xn, begin, begin + length / 2 - 1);
		FFT(xn, begin + length / 2, end);
	}

	Complex temp;
	for (int i = 0; i <= length / 2 - 1; i++)
	{
		temp = xn[begin + i];
		xn[begin + i + length / 2] = xn[begin + i + length / 2] *
			Complex(cos(2 * i * PI / length), -sin(2 * i * PI / length));//旋转因子
		xn[begin + i] = temp + xn[begin + i + length / 2];//基本蝶形计算
		xn[begin + i + length / 2] = temp - xn[begin + i + length / 2];
	}
	return xn;
}

vector<Complex> myFFT(vector<Complex> &xn)
{
	FFT(xn, 0, xn.size() - 1);
	return xn;
}

int main(void)
{
	vector<Complex> list = { 1, 1, 1, 1, 1, 1, 1, 1 };
	myFFT(list);

	for (auto it = list.begin(); it < list.end(); it++)//debug
	{
		cout << *it;
	}
	cout << endl;

	return 0;
}

运算结果:

1
(8.000+0.j)(0+0.j)(0+0.j)(0+0.j)(0+0.j)(0+0.j)(0+0.j)(0+0.j)

另外我在网上看到一个原生C接口的FFT函数库FFTW(the Fastest Fourier Transform in the West), 据说是世界上FFT的最快实现. Tutorial在这里.

下面是一个使用示例. 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#include "fftw3.h"
#include <iostream>
#pragma comment(lib, "libfftw3-3.lib") // double版本
// #pragma comment(lib, "libfftw3f-3.lib")// float版本
// #pragma comment(lib, "libfftw3l-3.lib")// long double版本

#define N 8

using namespace std;

int main(void)
{
	int i;
	fftw_complex *in, *out;
	fftw_plan p;
	in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
	out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);

	for (i = 0; i < N; i++) //输入序列
	{
		in[i][0] = 1;//实部
		in[i][1] = 0;//虚部
	}

	p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
	fftw_execute(p);

	for (i = 0; i < N; i++)
	{
		cout << "(" << in[i][0] << ", " << in[i][1] << "j)";
	}

	cout << endl;

	for (i = 0; i < N; i++)
	{
		cout << "(" << out[i][0] << ", " << out[i][1] << "j)";
	}

	cout << endl;

	fftw_destroy_plan(p);
	fftw_cleanup();
	fftw_free(in);
	fftw_free(out);

	return 0;
}

运算结果:

1
2
(1, 0j)(1, 0j)(1, 0j)(1, 0j)(1, 0j)(1, 0j)(1, 0j)(1, 0j)
(8, 0j)(0, 0j)(0, 0j)(0, 0j)(0, 0j)(0, 0j)(0, 0j)(0, 0j)

结果与我实现的FFT一样, 不过效率…

使用FFTW除需要先配置一下环境外, 还需要知道几个点:

  • 需添加#include "fftw3.h", 这不用说了吧.

  • 最好使用FFTW中实现的fftw_complex复数变量, 当如果添加了#include <complex>, 将会使用<complex>文件中定义的复数变量.

  • fftw_complex复数变量底层实现默认是个数组double[2], 依次存放了实虚部.

  • 另一个比较重要的是fftw_plan变量, 他就像是个机器, 把输入输出序列、序列长度以及其他一些参数传入, 使用fftw_execute(p)函数启动FFT运算.

  • 最后不能忘了清理资源, 务必使用fftw_free函数释放资源.

FFTW除了能够处理一维序列的DFT, 还可以处理二维序列及其他功能, 参见Tutorial.
.
.
.
.
.
.
.
.
.
終わり