seo网站编辑专业,长沙景点推荐,潍坊网站seo,服务器租用泰海传统的计算快速傅里叶变换的Cooley-Tukey算法效率极高#xff0c;因其主要由蝶形运算构成#xff0c;所以代码形式也非常简单#xff0c;只是需要将输入或者输出按照位反转的方式重新排序。这个重新排序的步骤并不是必须的。Clive Temperton于1991年在Self-Sorting In-Place…传统的计算快速傅里叶变换的Cooley-Tukey算法效率极高因其主要由蝶形运算构成所以代码形式也非常简单只是需要将输入或者输出按照位反转的方式重新排序。这个重新排序的步骤并不是必须的。Clive Temperton于1991年在Self-Sorting In-Place Fast Fourier Transforms一文中给出了适用于混合基数的原地FFT算法不需要对输入或输出重新排序。本文将介绍这种算法的原理并给出基数2Radix-2情况下的具体构造和C实现。作为FFT算法研究成果的集大成者FFTW已应用了这种算法。FFT7071专栏fourier.v.ariant.cn预备内存中的矩阵转置设x[t]为一个长度为M*N的向量t也可表示为Mab。以M5,N3为例x[0…14]{101, 102, 103,…, 115}如果将x视作一个5行3列的矩阵那么a列b行的矩阵元素即是x[Mab] a 0, 1, 2
b0 101 106 111
b1 102 107 112
b2 103 108 113
b3 104 109 114
b4 105 110 115将这个矩阵转置不难发现转置后的y[0…14]{101, 106, 111, 102,…, 110, 115} a 0, 1, 2, 3, 4
b0 101 102 103 104 105
b1 106 107 108 109 110
b2 111 112 113 114 115y与x的关系是y[Nba]x[Mab]。这个转置变换也可以用一个置换矩阵P表示yPx。展开DFT变换式记长度为 的信号为 为 的离散傅里叶变换并且设展开DFT变换得到以下表达式其中 。利用 用y代换x并继续展开单位根的幂其中 上式的求和等价于其中大括号内的求和相当于将y中的元素从地址0开始每相邻的N个为一组总共M个长度为N的DFT。如果将y看作M行N列的矩阵这是对每一列的变换变换的结果依次乘以 这时剩下的一个求和相当于对y的每一行单独的DFT。至此长度为MN的变换分解为了长度为M和N的两遍短变换。如果上式中不将第一遍DFT的结果乘以 结果将是M*N的2维DFT。需要注意的是变换y可以使上式的两遍DFT均在原地进行如果变换的是x为保持变换结果的顺序正确必须以转置的形式写回第一遍短变换的结果。题外话将中间结果写到另一处存储区x并且以x为输入做第二遍变换结果写回x如此往复可以解决变换无法在原地进行的问题这即是Stockham FFT算法。但是如此一来FFT需要额外的等于x长度的内存即需要额外O(N)的空间。因为置换群中的元素均可分解为2置换也就是对换的乘积置换矩阵也可做如此分解将转置操作变换成一系列对换从而可在原地进行仅需要O(1)额外空间。然而转置只能分解为数量巨大的对换这种操作的效率不高。这也预示着充分利用内存中数据的对换可以保持O(1)的额外空间需求同时使FFT在原地进行且顺序正确。展开DFT变换矩阵和素因子算法运用上文得出的分解到DFT的矩阵形式实际变换的是y也就是转置的x第一遍DFT作于相邻的N个元素每列将结果逐个乘以一个旋转因子再变换间隔为N的每组元素每行这个过程对应DFT矩阵的一种分解也就是素因子分解FFT算法的基本构造其中W为下标相应长度的DFT矩阵P为x[Mab]转置到y[Nba]的置换矩阵I为M或N阶的单位矩阵D为M*N阶对角矩阵第bNd行对角线上的值为exp(2πibd/(MN))。⊗是矩阵的Kronecker积定义如下例如Kronecker积满足结合律满足“混合乘积”性质以 为例继续分解 运用“混合乘积”的性质将上式拆分为矩阵积对于长度为 的DFT矩阵分解设可以得到这种分解正是时间抽取DIT基数2Radix-2的Cooley-Tukey算法下文中只考虑此种FFT频率抽取DIF在附录中讨论。已知T对应Cooley-Tukey算法每次迭代在整个输入向量上的所有蝶形运算上式中的 为2行8列到8行2列的矩阵转置作用是将x[2ba]的值变换到x[8ab]其中 。观察8ab和2ba的二进制表示2^3 2^2 2^1 2^0
[b] [b] [b] [a] 2ba
[a] [b] [b] [b] 8ab可知 的作用是将 的地址t二进制位向右环移1位。 的作用是保持t的最高1位不变t的余下三位向右环移1位因此经过所有的 变换2^3 2^2 2^1 2^0
[k3] [k2] [k1] [k0]
[k0] [k3] [k2] [k1] - P16
[k0] [k1] [k3] [k2] - P8
[k0] [k1] [k2] [k3] - P4很明显T之前所有 矩阵的乘积是输入数据翻转对应地址二进制位的置换矩阵。对换蝶形运算和对换设DFT的长度为N则x[t]中的t在二进制下有N位定义 为对换 和 的置换矩阵 由 对换二进制位中低位i和高位N-i-1得到。可以用一系列Q的乘积取代P的乘积 的效果同样完全反转二进制位2^3 2^2 2^1 2^0
[k3] [k2] [k1] [k0]
[k0] [k2] [k1] [k3] - Q03
[k0] [k1] [k2] [k3] - Q12在N2M或2M1的情况下转置P无法简单地在原地计算而Q仅包含数量较少的对换因此可以在原地完成。目前为止以T和Q组成的DFT矩阵W分解仍然表示先重排数据再开始蝶形运算的迭代如果将T和Q结合起来就能省略重排数据的操作但是这要求T和Q可交换。为了证明这一点首先需要求出Q的表达式。观察 翻转二进制最高位和最低位的操作0000 0000
0001 - 1000
0010 0010
0011 - 1010
0100 0100
0101 - 1100
0110 0110
0111 - 1110
1000 - 0001
1001 1001
1010 - 0011
1011 1011
1100 - 0101
1101 1101
1110 - 0111
1111 1111可以发现 的作用是将前一半数中的奇数 和 对换。因此这里 是 阶置换矩阵对换 和 由 交换二进制的最高位和最低位得到。在为二进制数添加前缀的操作下 的作用是不变的00000 00000 10000 10000
00001 - 01000 10001 - 11000
00010 00010 10010 10010
00011 - 01010 10011 - 11010
00100 00100 10100 10100
00101 - 01100 10101 - 11100
00110 00110 10110 10110
00111 - 01110 10111 - 11110
01000 - 00001 11000 - 10001
01001 01001 11001 11001
01010 - 00011 11010 - 10011
01011 01011 11011 11011
01100 - 00101 11100 - 10101
01101 01101 11101 11101
01110 - 00111 11110 - 10111
01111 01111 11111 11111因此对于N位二进制数为二进制数添加后缀的操作使 作用于全部后缀可得出设 现在可以将 展开因此对于所有的 和 可交换。这使 可以写为已知一次蝶形运算的迭代 的作用是将两个长度为 的DFT结果作为奇偶两部分合并为长度为 的DFT结果这样的奇偶对共有 个。令 中单个蝶形运算的偶、奇两个输入分别是 和 会将 与 对换。取 则 将与 对换。在 的情况下 与 分别是另一蝶形运算的偶、奇输入。可见 中输入数据地址相差 的两个蝶形运算是成对的第一个蝶形运算的奇数项输入与第二个蝶形运算的偶数项输入对换。如下图所示下图中画出来N16的基数2变换输入和输出的顺序均是正确的图中用颜色标出了某些蝶形运算使蝶形运算的配对更清晰。注意其中成对蝶形运算的4个输入输出均在原地并且与传统Cooley-Tukey算法相比没有计算量的差别。下图是作为对比的传统Cooley-Tukey算法。附录本文算法的C实现/* copyright 2020, github.com/zhxxch, all rights reserved. */
#include complex
#include iterator
#include cmath
#include cassert
#include cstddef
templatetypename iter_t
#if __cplusplus 201703L
requires std::random_access_iteratoriter_t
#endifinline void fft_in_place(iter_t arr_begin,iter_t arr_end, const bool is_inverse) {using cplx_t typename std::iterator_traitsiter_t::value_type;using real_t typename cplx_t::value_type;const size_t length std::distance(arr_begin, arr_end);assert(requires length pow(2,n) (length (length - 1)) 0);constexpr real_t pi 3.141592653589793238462643383;size_t sub_ft_size 1;size_t num_sub_ft length / sub_ft_size;size_t num_sub_ft_pair num_sub_ft / 2;for(; sub_ft_size num_sub_ft_pair;sub_ft_size * 2, num_sub_ft / 2,num_sub_ft_pair / 2) {for(size_t coupled_group_pos 0;coupled_group_pos length;coupled_group_pos 2 * num_sub_ft_pair) {for(size_t sub_ft_pos coupled_group_pos;sub_ft_pos coupled_group_pos num_sub_ft_pair;sub_ft_pos 2 * sub_ft_size) {for(size_t i sub_ft_pos, nth_pow 0;i sub_ft_pos sub_ft_size;i, nth_pow num_sub_ft_pair) {const cplx_t W exp((is_inverse ? 1. : -1.)* cplx_t(0, 2 * nth_pow * pi)/ (real_t)length);const cplx_t parit00 arr_begin[i];const cplx_t parit01 arr_begin[num_sub_ft_pair i]* W;const cplx_t parit10 arr_begin[i sub_ft_size];const cplx_t parit11 arr_begin[num_sub_ft_pair i sub_ft_size]* W;arr_begin[i] parit00 parit01;arr_begin[i sub_ft_size] parit00 - parit01;arr_begin[num_sub_ft_pair i] parit10 parit11;arr_begin[num_sub_ft_pair i sub_ft_size] parit10 - parit11;}}}}for(; sub_ft_size length; sub_ft_size * 2,num_sub_ft / 2, num_sub_ft_pair / 2) {for(size_t sub_ft_pos 0; sub_ft_pos length;sub_ft_pos 2 * sub_ft_size) {for(size_t i sub_ft_pos, nth_pow 0;i sub_ft_pos sub_ft_size;i, nth_pow num_sub_ft_pair) {const cplx_t parit1 arr_begin[i sub_ft_size]* exp((is_inverse ? 1. : -1.)* cplx_t(0, 2 * nth_pow * pi)/ (real_t)length);const cplx_t parit0 arr_begin[i];arr_begin[i] parit0 parit1;arr_begin[i sub_ft_size] parit0 - parit1;}}}
}
使用方法-FFTfft_in_place(arr.begin(), arr.end(), 0);
使用方法-IFFTfft_in_place(arr.begin(), arr.end(), 1);
arr的长度必须是2的幂。附录时间抽取与频率抽取以 为例频率抽取的FFT算法中换为使用Q表达的形式则为因此Q作用于蝶形运算的输出。
