Pffft:a pretty fast FFT。
1.官网下载
https://bitbucket.org/jpommier/pffft/src/master/
pffft.h和pffft.c就是pffft的全部代码,fftpack.h和fftpack.c是fftpack库,test_pffft.c是对比测试pffft和fftpack的。结果当然是pffft比fftpack运行速度快,但是在不使用SIMD时,二者运行速度几乎是一样的;但开启SIMD后,pffft的速度明显就快很多了。
2.VS2015的编译
Pffft默认是开启SIMD的,但是需要对VS做一下配置。
- 优化
- 预处理器
- 代码生成
在Vs下,使用SIMD,只有Release模式下能编译过,Debug模式编译不成功。如果不使用SIMD,需要预处理器定义中添加PFFFT_SIMD_DISABLE。
3.复数例子
typedef float pffft_scalar;
int complex_test(int N, int cplx, int useOrdered)
{
pffft_scalar *X = pffft_aligned_malloc((unsigned)N * 2 * sizeof(pffft_scalar));
pffft_scalar *Y = pffft_aligned_malloc((unsigned)N * 2 * sizeof(pffft_scalar));
pffft_scalar *R = pffft_aligned_malloc((unsigned)N * 2 * sizeof(pffft_scalar));
pffft_scalar *Z = pffft_aligned_malloc((unsigned)N * 2 * sizeof(pffft_scalar));
pffft_scalar *W = pffft_aligned_malloc((unsigned)N * 2 * sizeof(pffft_scalar));
//assert(pffft_is_power_of_two(N));
PFFFT_Setup *s = pffft_new_setup(N, PFFFT_COMPLEX);
assert(s);
if (!s) {
printf("Error setting up PFFFT!\n");
return 1;
}
for (int k = 0; k < 2 * N; k += 2)
{
X[k] = 1.0f* (k / 2) + 1;
X[k + 1] = X[k];
}
pffft_transform_ordered(s, X, Y, W, PFFFT_FORWARD);//fft
pffft_transform_ordered(s, Y, Z, W, PFFFT_BACKWARD);//ifft
//Z是ifft的结果,需要除以N
pffft_destroy_setup(s);
pffft_aligned_free(X);
pffft_aligned_free(Y);
pffft_aligned_free(Z);
pffft_aligned_free(R);
pffft_aligned_free(W);
return 1;
}
4.实数例子
int real_test(int N, int cplx, int useOrdered)
{
pffft_scalar *X = pffft_aligned_malloc((unsigned)N * sizeof(pffft_scalar));
pffft_scalar *Y = pffft_aligned_malloc((unsigned)N * sizeof(pffft_scalar));
pffft_scalar *R = pffft_aligned_malloc((unsigned)N * sizeof(pffft_scalar));
pffft_scalar *Z = pffft_aligned_malloc((unsigned)N * sizeof(pffft_scalar));
pffft_scalar *W = pffft_aligned_malloc((unsigned)N * sizeof(pffft_scalar));
pffft_scalar *Y2 = pffft_aligned_malloc((unsigned)N *2* sizeof(pffft_scalar));
//assert(pffft_is_power_of_two(N));
PFFFT_Setup *s = pffft_new_setup(N, PFFFT_REAL);
assert(s);
if (!s) {
printf("Error setting up PFFFT!\n");
return 1;
}
for (int k = 0; k < N; k++){
X[k] = 1.0f* k;
}
double t0, t1, flops;
t0 = uclock_sec();
pffft_transform_ordered(s, X, Y, W, PFFFT_FORWARD);
for (int i = 0; i < N; i++) {
printf("%d: %f \n", i, Y[i]);
}
printf("==================\nfft:\n");
{
for (int i = 2; i <= N; i++){
Y2[i] = Y[i];
}
Y2[0] = Y[0];
Y2[1] = 0;
Y2[N] = Y[1];
Y2[N + 1] = 0;
for (int i = N + 2; i <= 2 * N - 1; i++)
{
if (i % 2)
Y2[i] = -Y2[2 * N - i + 2];
else
Y2[i] = Y2[2 * N - i];
}
}
for (int i = 0; i < 2 * N; i += 2)
printf("%d %f %f\n", i / 2, Y2[i], Y2[i + 1]);
printf("======================\nifft:\n");
//printf("=========================\n");
//t0 = uclock_sec();
//if (useOrdered)
pffft_transform_ordered(s, Y, Z, W, PFFFT_BACKWARD);
t1 = uclock_sec();
//else
//pffft_transform(s, R, Z, W, PFFFT_BACKWARD);
//for (int i = 0; i < N; i++) {
// printf("%d: %f %f\n", i, Y[2 * i], Y[2 * i + 1]);
//}
//printf("=========================\n");
for (int i = 0; i < N; i++) {
printf("%d: %f \n", i, Z[i]/N);
}
printf("pffft_transform_ordered PFFFT_BACKWARD %d %f \n", N, t1 - t0);
pffft_destroy_setup(s);
pffft_aligned_free(X);
pffft_aligned_free(Y);
pffft_aligned_free(Z);
pffft_aligned_free(R);
pffft_aligned_free(W);
pffft_aligned_free(Y2);
return 1;
}
5.分析
- 与matlab结果比对一致,运算速度还可以,比自己写的fft快了很多。
- 使用VS,Debug模式不支持SIMD。
- Pffft的精度是float,不是double。Github上有double版本的pffft,但是不支持SIMD。