fftw3笔记_fftw_plan_dft_1d

微信扫一扫,关注公众号“音频算法与工程实践”

Linux:http://www.fftw.org/download.html
Windows:http://www.fftw.org/install/windows.html
old version:ftp://ftp.fftw.org/pub/fftw/
---------------------------------------------app代码-----------------------------------
原文参考fftw3.pdf的"2 Tutorial", "4.1.2 Precision"
1.FFTW有三个版本的数据类型float, double, long double。默认是double版本
2.都使用同样的头文件fftw3.h
3.所有以小写fftw_开头的比如函数和类型替换为fftwf_或者fftwl_, fftw_complex 改为 fftwf_complex
以大写FFTW_开头的保持不变。
4.函数参数中double替换成float或者long double
5.fftw_complex, which is by default a double[2] composed of the real (in[i][0]) and imaginary (in[i][1])
6.一维复数据的DFT
fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags);
如果in和out指针相同为原位运算，否则为非原位运算。输入in有n个数据，输出out有n个数据。
sign可以为正变换 FFTW_FORWARD(-1)，也可以为逆变换 FFTW_BACKWARD(+1)，实际上就是变换公式中指数项的符号。
一维实数据的DFT
fftw_plan fftw_plan_dft_r2c_1d(int n, double *in, fftw_complex *out, unsigned flags);
fftw_plan fftw_plan_dft_c2r_1d(int n, fftw_complex *in, double *out, unsigned flags);
r2c版本：实输入数据，复Hermitian输出，正变换。输入in有n个数据，输出out有floor(n/2)＋1个数据。
c2r版本：复Hermitian输入数据，实输出数据，逆变换。
一个fftw_plan只能对长度固定的in/out进行变换，但可以在变换后改变in的内容（大小不变）以用同一个方案执行新的变换。
实数据DFT中，首个变换数据为直流分量，为实数；如果n为偶数，第n/2个变换数据也为实数；所以可以把这两个数据组合在一起，即将第n/2个变换数据作为第0个变换数据的虚部。有些FFT的实现就是这么做的，FFTW没有这么做。
7.需注意FFTW的逆变换没有除以N，即数据正变换再反变换后是原始数据的N倍。
---------------------------------------------编译lib------------------
./configure --enable-float --disable-fortran
make
make install
"make install" command typically requires root privileges.

fftw-3.3.2.tar.g:VS2017编译
unzip fftw-3.3.2.tar.gz to D:/fftw-3.3.2
unzip fftw-3.3-libs-visual-studio-2010.zip to D:/fftw-3.3.2/fftw-3.3-libs
open fftw-3.3-libs.sln, SDK choose 10.0; Platform choose VS2017; Include Directories add $(ProjectDir)\..\..;
delete \libbench2\aligned-main.c
Rebuild Solution.
D:\fftw-3.3.2\fftw-3.3-libs\x64\Debug

fftw-3.3.5-dll32.zip:VS2017生成lib
cd C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Tools\MSVC\14.13.26128\bin\Hostx64\x64\ (Hostx64\x86\lib.exe文件不全)
copy libfftw3f-3.def ./
lib /machine:x86 /def:libfftw3f-3.def
官网和"README-WINDOWS"都是使用"lib /def:libfftw3f-3.def"，应该与VC版本有关。
---------------------------------------------编译app----------------------------------------
https://blog.csdn.net/wangwei_cq/article/details/8187576 提到链接.a文件编译app，但官方下载没有该文件，只能先将def转换成lib，再链接dll和lib两个文件。
Build options - test - Linker settings - Add 选择 bin\Debug\libsimple.a - 确定，编译即可成功。
CodeBlocks编译app需要添加fftw3.h和libfftw3l-3.lib，运行时还需要dll文件。默认是x86。
VS2017编译app方法类似
1.将头文件装到#include目录下
2.将库文件装到#lib目录下
3.将动态链接库装到系统Windows32目录下

window下查看dll, lib文件是32位还是64位命令:
cd C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Tools\MSVC\14.13.26128\bin\Hostx64\x64

dumpbin.exe /headers libfftw3f-3.lib

------------------------------------------------demo----------------------------

#define NN 8

#include "fftw3.h"
//fftw-3.3.5-dll32
// https://blog.csdn.net/tuzanbo/article/details/75309647
void fftw3_test()
{
    int len = NN;
    float *in = NULL;
    int i;
    // 如果要使用float版本,需先引用float版本的lib库,然后在fftw后面加上f后缀即可.
    fftwf_complex *out = NULL; // fftwf_complex --> 即为float版本
    fftwf_plan p;
    fftwf_plan pifft;
    in = (float *) fftwf_malloc(sizeof(float) * len);
    out = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex) * len);
    //float dx = 1.0 / len;


    for(i = 0; i < len; i++)
    {
        in[i] = i+1;
        printf("%f ", in[i]);
    }
    printf("\n\n");


    // 傅里叶变换
    p = fftwf_plan_dft_r2c_1d(len, in, out, FFTW_ESTIMATE);
    fftwf_execute(p);
    for(i = 0; i < len; i++)
{
        printf("%f%+fi\n", out[i][0], out[i][1]);
    }
    printf("\n");


    // 输出幅度谱
    for(i = 0; i < len; i++)
{
        float len = sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]);
        printf("%.2f ", len);
    }
    printf("\n");


    pifft = fftwf_plan_dft_c2r_1d(len, out, in, FFTW_ESTIMATE);
    fftwf_execute(pifft);
    for(i = 0; i < len; i++)
{
        printf("%f %f\n", in[i], in[i]/len);
    }
    printf("\n");


    // 释放资源
    fftwf_destroy_plan(p);
    fftwf_destroy_plan(pifft);
    fftwf_free(in);
    fftwf_free(out);


    printf("-------------- end ---------------\n");
    return ;
}


void fftw3_dft_test()
{
    int len = NN;
    fftwf_complex *in = NULL;
    int i;
    // 如果要使用float版本,需先引用float版本的lib库,然后在fftw后面加上f后缀即可.
    fftwf_complex *out = NULL; // fftwf_complex --> 即为float版本
    fftwf_plan p;
    fftwf_plan pifft;
    in = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex) * len);
    out = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex) * len);
    //float dx = 1.0 / len;


    for(i = 0; i < len; i++)
    {
        in[i][0] = i+1;
        in[i][1] = 0;
        printf("%f%+fi\n", in[i][0], in[i][1]);
    }
    printf("\n\n");


    // 傅里叶变换
    p = fftwf_plan_dft_1d(len, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftwf_execute(p);
    for(i = 0; i < len; i++)
{
        printf("%f%+fi\n", out[i][0], out[i][1]);
    }
    printf("\n");


    // 输出幅度谱
    for(i = 0; i < len; i++)
{
        float len = sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]);
        printf("%.2f ", len);
    }
    printf("\n");


    pifft = fftwf_plan_dft_1d(len, out, in, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftwf_execute(pifft);
    for(i = 0; i < len; i++)
{
        printf("%f %f\n", in[i][0], in[i][0]/len);
    }
    printf("\n");


    // 释放资源
    fftwf_destroy_plan(p);
    fftwf_destroy_plan(pifft);
    fftwf_free(in);
    fftwf_free(out);


    printf("-------------- end ---------------\n");
    return ;
}

更多音频文章，请关注微信公众号“音频算法与工程实践”。