Pffft库的学习

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.分析

1. 与matlab结果比对一致,运算速度还可以,比自己写的fft快了很多。
2. 使用VS,Debug模式不支持SIMD。
3. Pffft的精度是float,不是double。Github上有double版本的pffft，但是不支持SIMD。