多维复数据的DFT同一维复数据的DFT用法类似，数组in/out为行优先方式存储。fftw_plan_dft是一个通用的复DFT函数，可以执行任意维复DFT。比如对于图像(2维数据)，其变换fftw_plan_dft_2d(height,width,85),因为原始图像数据为height×width的矩阵，即第一维(n0)为行数
height。
3.
一维实数据的DFT
在许多应用中，输入数据in是纯实数，此时DFT输出满足Hermitian冗余：out[i]
是out[n-i]的共轭。利用这个性质可以从速度和存储两方面改进算法。
实数的fftw调用过程与复数的基本相同，差别只在in/out是double型并且创建方案是用以下函数：
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输出，正变换。
c2r版本：复Hermitian输入数据，实输出，逆变换。
对于单精度或长双精度版本，double型应替换为float型或long
double型。
n：逻辑长度，不必为数组的物理长度。由于实数据的DFT具有Hermitian对称性，所以只需要计算n/2+1(除法向下取整)个复数输出就可以了。比如对于r2c，输入in有n个数据，输出out有floor(n
/2)＋1个数据。对于原位运算，in和out为同一数组(out须强制类型转换)，必须足够大以容纳所有数据，因而长度为2*(n/2+1)，数组前n个之后的元素未用，是给输出填充的。
c2r逆变换在任何情况下(不管是否为原位运算)都破坏输入数组in，如果有必要可以通过在flags中加入FFTW_PRESERVE_INPUT来阻止，但这会损失一些性能。这个标志位目前在多维实数DFT中已不被支持。
比如对于n=4，in=[1 2 3 4]，out=[10 -2+2i -2
-2-2i]，out具有共轭对称性，out的实际内存为10 0 -2 2 -2 0，共3个复数数据。对于n=5，in=[1 2 3 4
5]，out=[15 -2.5+3.44i -2.5+0.81i -2.5-0.81i -2.5-3.44i]
，out的实际内存为15 0 -2.5 3.44 -2.5 0.81，共3个复数数据。
实数据DFT中，第0个变换数据为DC分量，为实数；如果n为偶数，由
Nyquist采样定理，第n/2个变换数据也为实数；所以可以把这两个数据组合在一起，即将第n/2个变换数据作为第0个变换数据的虚部，这样一来输入数组就和输出数组的个数相等(n=n/2*2)。一些FFT的实现就是这么做的，但FFTW没有这么做，因为它并不能很好地推广到多维DFT中，而且存储空间的节省也是非常小以至于可以忽略不计。
一个一维r2c和c2r
DFT的替代接口可以在r2r接口中找到，它利用了半复数输出类型(即实部和虚部分开放在不同的区域里)，使输出数组具有和输入数组同样的长度和类型。该接口尽管在多维变换中用处不大，但有时可能会有一些性能的提升。
4.多维实数的DFT
fftw_planfftw_plan_dft_r2c_2d( int n0, int n1, double *in, fftw_complex *out, unsigned flags);
fftw_plan fftw_plan_dft_r2c_3d( int n0, int n1, int n2, double *in, fftw_complex *out, unsigned flags);
fftw_plan fftw_plan_dft_r2c( int rank, const int *n, double *in, fftw_complex *out, unsigned flags);
这是r2c接口(正变换)，c2r接口(逆变换)只是简单的将输入/输出类型交换一下。用法大致同1d情况一样，需要特别注意的是复数据的存放方式。对于n0×n1×n2×…×nd-1的实输入数组(行优先)，经过r2c正变换后，输出为一个n0×n1×n2×…×(nd-1/2+1)的
fftw_complex
型复数数组(行优先)，其中除法向下取整。由于复数数据的总长度要大于实数据，所以如果需要进行原位运算则输入实数组必须扩展以能够容纳所有复数据，即实数组的最后一维必须包含2(floor(nd-1/2)+1)个double元素。比如，对于一个2维实正DFT，输入数组为n0×n1大小，输出复数组大小为n0×floor(n1/2+1)(由对称性)，其需要的内存大于实输入数组。所以对于原位运算，输入数组要扩展到n0×2floor(n1/2+1)。如果n1为偶数，扩展为n0×(n1＋2)；如果n1为奇数，扩展为n0×(n1＋1)；这些扩展的内存不需要赋初值，它们只用来存放输出数据。
比如对于3×3的实正DFT，in=[0
2 4; 6 1 3; 5 7 4]，out=[32 0.5+0.86i 0.5-0.86i; -7+5.2i -1-1.73i
-8.5-6.06i; -7-5.2i -8.5+6.06i
-1+1.73i]；out的实际内存为32,0,0.5,0.86,-7,5.2,-1,-1.73,-7,-5.19,-8.5,6.06；即为
3×2的复数组，换算成double类型为3×4，所以如果要进行原位运算，in数组大小必须为3×4，即最后一维(每行)扩展一个double元素。另外，如果用这个out数组进行3×3的c2r逆变换，将得到实数据[0
18 36; 54 9 27; 45 63
36]，即原始数据的9倍，这是因为FFTW的逆变换没有归一化。
5.
更多实数的DFT
通过一个统一的r2r(real-to-real)接口，FFTW支持其它的一些变换类型，这些变换的输入和输出数组都为double型，且大小相同。这些r2r变换可以分为3个类型：实输入，复Hermitian对称输出(以halfcomplex格式)的DFT；DCT/DST(偶/奇对称的实输入DFT，即离散余弦/正弦变换)；DHT(离散
Hartley变换)。
fftw_planfftw_plan_r2r_1d(int n, double *in, double *out, fftw_r2r_kind kind, unsigned flags);
fftw_plan fftw_plan_r2r_2d(int n0, int n1, double *in, double *out, fftw_r2r_kindkind0, fftw_r2r_kind kind1, unsigned flags);
fftw_plan fftw_plan_r2r_3d(int n0, int n1, int n2, double *in, double *out, fftw_r2r_kind kind0, fftw_r2r_kind kind1, fftw_r2r_kind kind2, unsigned flags);
fftw_plan fftw_plan_r2r(int rank, const int *n, double *in, double *out, const fftw_r2r_kind *kind, unsigned flags);
这里n为数组的物理尺寸。对于多维变换，数组按行优先方式存储(与C++标准相同，与Fortran不同)。由于DFT是可分离变换，所以2维/3维/多维的变换是在每个维度上分别进行变换得到的，每个维度都可指定一个kind参数，指定该维的r2r变换类型。
(1).
半复数格式(HalfComplex-format)
一种r2r的kind是FFTW_R2HC (r2hc)
对应r2c形式DFT，但是以半复数格式作为输出，有时这种格式更为快速和方便。它的逆变换是kind为FFTW_HC2R
(hc2r)。对于大小为n的1维实输入的DFT，输出n个double型实数，格式如下：
r0,
r1, r2, … , rn/2,
i(n+1)/2-1, … , i2,
i1
rk是第k个输出的实部，ik是第k个输出的虚部。(n+1)/2向下取整。对于一个半复数格式的数组hc[n]，第k个元素的实部为hc[k]，虚部为hc[n-k]，除了k==0或n/2(n为偶数)的情况，这两种情况下由于实输入的对称性，虚部为0，不存储。所以对于r2hc
正变换，输入输出数组大小都为n，hc2r
逆变换也是如此。除了数据格式的差异，r2hc和hc2r变换的结果与前述的1维r2c和c2r的变换结果是相同的。
对于多维比如2维变换，由可分离性，可以先对行变换，再对列变换，FFTW_R2HC方式行变换的结果是半复数行，如果采用FFTW_R2HC方式进行列变换，需要进行一些数据处理，r2r变换不会自动处理，需要手动进行，所以对于多维实数据变换，推荐使用普通的r2c/c2r接口。
(2). DCT/DST
离散余弦变换DCT是实偶对称数据的DFT(REDFT, Real-Even DFT),
离散正弦变换DST是实奇对称数据的DFT(RODFT, Real-Odd
DFT)。REDFTab和RODFTab中的a,b为数据移位标志，表示输入a和输出b是否被移位(1表示移位)，这些构成了DCT和DST的I－IV类，比较常用的为DCT-II，FFTW的r2r接口支持所有四种类型的变换。对于n个实数的输入数组in[
j=0, … , n-1]，
FFTW_REDFT00 (DCT-I): even around j=0 and even
around j=n-1.
FFTW_REDFT10 (DCT-II, the DCT): even around
j=-0.5 and even around j=n-0.5.
FFTW_REDFT01 (DCT-III, the IDCT): even around j=0
and odd around j=n.
FFTW_REDFT11 (DCT-IV): even around j=-0.5 and odd
around j=n-0.5.
FFTW_RODFT00 (DST-I): odd around j=-1 and odd
around j=n.
FFTW_RODFT10 (DST-II): odd around j=-0.5 and odd
around j=n-0.5.
FFTW_RODFT01 (DST-III): odd around j=-1 and even
around j=n-1.
FFTW_RODFT11 (DST-IV): odd around j=-0.5 and even
around j=n-0.5.
注意对称性只是逻辑意义上的，对物理输入数据没有任何限制。比如，对于n＝5的REDFT00
(DCT-I)，输入数据为abcde，它对应DFT输出为n=8的逻辑偶数组abcdedcb；对于n＝4的REDFT10
(DCT-II)，输入数据为abcd，它对应n=8的逻辑偶数组abcddcba。
所有这些变换都是可逆的。R*DFT00的逆变换是R*DFT00，R*DFT10的逆变换是R*DFT01(即通常所说的DCT和IDCT)，R*DFT11的逆变换是R*DFT11。如同实数和复数的DFT一样，这些变换的结果都没有归一化，所以正变换再逆变换后数据变为原始数据的N倍，N为逻辑的DFT大小。对于REDFT00，N=2(n-1)；对于RODFT00，N=2(n+1)；否则N=2n。
对于如何选择哪一种变换的问题，需要根据实际问题确定，不过有以下一些建议。R*DFT01和R*DFT10要比R*DFT11稍微快一些，尤其对于奇数长度数据；而R*DFT00有时要更慢一些，尤其对于奇数长度数据。
比如对于in=[1 2 3
4]，n=4，N=2n=8。Matlab下dct变换的结果为[5 -2.2304 0
-0.15851]；FFTW的结果为(FFTW_REDFT10)out=[20 -6.3086 0
-0.]，为matlab结果的√8(√N)倍；用out进行逆dct变换(FFTW_REDFT01)的结果为[8 16
24 32]，为原始数据的8(N)倍。
再比如对于in=[0 2 4; 6 1
3; 5
7 4]的二维DCT变换，n=3，N=2n=6。Matlab下dct2的变换结果为out_matlab=[10.667
0 0.4714; -4.0825 -2.5 1.4434; 0.4714 -2.5981
-3.1667]；FFTW的结果为(fftw_plan_r2r_2d(3, 3, in, out, FFTW_REDFT10,
FFTW_REDFT10, FFTW_ESTIMATE)out_fftw=[128 0 4; -34.641 -15 8.66; 4
-15.588
-19]，这与Matlab的结果同样是有差别的。将Matlab的结果变换到FFTW结果的步骤为：
1).
将out_matlab乘以√6×√6(即√N×√N)；
2).
再将上一步得到的out_matlab的第一行和第一列都乘以√2，因此第一个元素(即左上角的元素)要乘以2。
第一个是归一化的原因，第二个是FFTW为了将DCT变换与实偶对称FFT相对应的结果。这些对于DCT变换的应用都影响不大。
最后对out_fftw进行IDCT变换
(fftw_plan_r2r_2d(3, 3, in, out, FFTW_REDFT01, FFTW_REDFT01,
FFTW_ESTIMATE)，将得到[0 72 144; 216 36 108; 180 252
144]；它是原始数据in的36(N×N，N=6)倍。
6.
数据类型
FFTW有三个版本的数据类型：double、float和long
double，使用方法如下：
链接对应的库(比如libfftw3-3、libfftw3f-3、或ibfftw3l-3)
包含同样的头文件fftw3.h
将所有以小写”fftw_”开头的名字替换为”fftwf_”(float版本)或”fftwl_”(long
double版本)。
所有以大写”FFTW_”开头的名字不变
将函数参数中的double替换为float或long
double
最后，虽然long
double是C99的标准，但你的编译器可能根本不支持该类型，或它并不能提供比double更高的精度。
fftw_malloc考虑了数据对齐，以便使用SIMD指令加速，所以最好不要用C函数malloc替代，而且不要将fftw_malloc、fftw_free和malloc、free、new、delete等混用。尽量使用fftw_malloc分配数组，而不是下面的固定数组，因为固定数组是在栈上分配的，而栈空间较小；还因为这种方式没有考虑数据对齐，不便应用SIMD指令。
还有2个封装函数 double*
fftw_alloc_real(size_t n) 和 fftw_complex* fftw_alloc_complex(size_t
n) 便于使用，等价于实数和复数形式的(double*) fftw_malloc(sizeof(double)*n) 和
(fftw_complex *)
fftw_malloc(sizeof(fftw_complex)*n)。
7.基本函数接口
FFTW的API分为三类：基本接口，高级接口和guru接口。这里只介绍基本接口。
1).
复数DFT
fftw_planfftw_plan_dft_1d(int n,
fftw_complex *in, fftw_complex *out,
int sign, unsigned flags);
fftw_planfftw_plan_dft_2d(int n0, int n1,
fftw_complex *in, fftw_complex *out,
int sign, unsigned flags);
fftw_planfftw_plan_dft_3d(int n0, int n1, int n2,
fftw_complex *in, fftw_complex *out,
int sign, unsigned flags);
fftw_planfftw_plan_dft(int rank, const int *n,
fftw_complex *in, fftw_complex *out,
int sign, unsigned flags);
2).
实数DFT
fftw_planfftw_plan_dft_r2c_1d(int n,
double *in, fftw_complex *out,
unsigned flags);
fftw_planfftw_plan_dft_r2c_2d(int n0, int n1,
double *in, fftw_complex *out,
unsigned flags);
fftw_planfftw_plan_dft_r2c_3d(int n0, int n1, int n2,
double *in, fftw_complex *out,
unsigned flags);
fftw_planfftw_plan_dft_r2c(int rank, const int *n,
double *in, fftw_complex *out,
unsigned flags);
以上是正变换，它们的逆变换形式是c2r，输入是复数fftw_complex型，输出是double型。
3).
实数－实数变换
fftw_planfftw_plan_r2r_1d(int n, double *in, double *out,
fftw_r2r_kind kind, unsigned flags);
fftw_planfftw_plan_r2r_2d(int n0, int n1, double *in, double *out,
fftw_r2r_kind kind0, fftw_r2r_kind kind1, unsigned flags);
fftw_planfftw_plan_r2r_3d(int n0, int n1, int n2,
double *in, double *out,
fftw_r2r_kind kind0,
fftw_r2r_kind kind1,
fftw_r2r_kind kind2,
unsigned flags);
fftw_planfftw_plan_r2r(int rank, const int *n, double *in, double *out, const fftw_r2r_kind *kind, unsigned flags);
4).
实数－实数变换类型
对于大小为n的下列变换，对应的逻辑DFT大小为N，N用来进行归一化。FFTW的变换没有归一化，正变换后再逆变换为原数据的N倍(不是n倍)，对于多维变换，为N的乘积(N0*N1*N285)。下表列出了变换类型及其对应的逻辑大小N及逆变换：
FFTW_R2HC computes a real-input DFT with output
in halfcomplex format, i.e. real and imaginary parts for a
transform of size n stored as:r0, r1,
r2, …, rn/2, i(n+1)/2-1, …,
i2, i1 (Logical N=n, inverse is
FFTW_HC2R.)
FFTW_HC2R computes the reverse of FFTW_R2HC,
above. (Logical N=n, inverse is FFTW_R2HC.)
FFTW_DHT computes a discrete Hartley transform.
(Logical N=n, inverse is FFTW_DHT.)
FFTW_REDFT00 computes an REDFT00 transform, i.e.
a DCT-I. (Logical N=2*(n-1), inverse is
FFTW_REDFT00.)
FFTW_REDFT10 computes an REDFT10 transform, i.e.
a DCT-II (sometimes called the DCT). (Logical N=2*n, inverse is
FFTW_REDFT01.)
FFTW_REDFT01 computes an REDFT01 transform, i.e.
a DCT-III (sometimes called the IDCT, being the inverse of DCT-II).
(Logical N=2*n, inverse is FFTW_REDFT=10.)
FFTW_REDFT11 computes an REDFT11 transform, i.e.
a DCT-IV. (Logical N=2*n, inverse is
FFTW_REDFT11.)
FFTW_RODFT00 computes an RODFT00 transform, i.e.
a DST-I. (Logical N=2*(n+1), inverse is
FFTW_RODFT00.)
FFTW_RODFT10 computes an RODFT10 transform, i.e.
a DST-II. (Logical N=2*n, inverse is
FFTW_RODFT01.)
FFTW_RODFT01 computes an RODFT01 transform, i.e.
a DST-III. (Logical N=2*n, inverse is
FFTW_RODFT=10.)
FFTW_RODFT11 computes an RODFT11 transform, i.e.
a DST-IV. (Logical N=2*n, inverse is
FFTW_RODFT11.)
8.
用同一个fftw_plan执行多个数据的变换
前面说过可以利用同一个fftw_plan通过对输入数据赋不同值来实现不同的变换，实际上还可以利用同一个fftw_plan实现对不同输入输出数据的变换，也就是说可以有多个输入输出数据数组，各自进行变换，互不影响。当然这要满足一定的条件：
输入/输出数据大小相等。
变换类型、是否原位运算不变。
对split数组(指实虚部分开)，实部和虚部的分割方式与方案创建时相同。
数组的对齐方式相同。如果都是用fftw_malloc分配的则此项条件满足，除非使用FFTW_UNALIGNED标志。
如果想对新数组，比如大小相等的一批数组执行变换，可以使用以下接口：
void fftw_execute_dft(
const fftw_plan p,
fftw_complex *in, fftw_complex *out);
void fftw_execute_split_dft(
const fftw_plan p,
double *ri, double *ii, double *ro, double *io);
void fftw_execute_dft_r2c(
const fftw_plan p,
double *in, fftw_complex *out);
void fftw_execute_split_dft_r2c(
const fftw_plan p,
double *in, double *ro, double *io);
void fftw_execute_dft_c2r(
const fftw_plan p,
fftw_complex *in, double *out);
void fftw_execute_split_dft_c2r(
const fftw_plan p,
double *ri, double *ii, double *out);
void fftw_execute_r2r(
const fftw_plan p,
double *in, double *out);
这些函数的执行不会修改原始plan，并且可以和fftw_execute混合使用。
9.
FFTW变换公式
一维复DFT正变换
一维复DFT逆变换
一维实DFT正变换(Yk = Yn-k*)
一维实DFT逆变换
REDFT00 (DCT-I)
REDFT10 (DCT-II)
REDFT01 (DCT-III)
REDFT11 (DCT-IV)
RODFT00 (DST-I)
RODFT10 (DST-II)
RODFT01 (DST-III)
RODFT11 (DST-IV)
1d Discrete Hartley Transforms (DHTs)
发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/225496.html原文链接：https://javaforall.net
