GSL中的快速傅里叶变换（FFTS）

快速傅里叶变换（FFTS）

    本章描述了执行快速傅里叶变换(FFTs)的函数。库中包含了radix-2例程(长度为2的幂)和混合基数例程(长度不限)。为了提高效率，对于真实数据和复杂数据，有不同版本的例程。混合基数例程是对Paul Swarztrauber的FFTPACK库的重新实现。FFTPACK的Fortran代码在Netlib上是可用的(FFTPACK也包括一些用于正弦和余弦变换的例程，但这些例程目前在GSL中不可用)。有关底层算法的详细信息和派生，请参阅文档“GSL FFT算法”(参见参考资料和进一步阅读)。

16.1 数学定义

    快速傅里叶变换是计算离散傅里叶变换(DFT)的有效算法，

    当函数在空间或时间上以离散间隔采样时，DFT通常是连续傅里叶变换的近似值。离散傅立叶变换的简单评估是矩阵向量乘积Wz。

一般的矩阵向量乘法对n个数据点进行O(n2)次运算。快速傅里叶变换算法使用分治策略将矩阵W分解成更小的子矩阵，这些子矩阵对应于长度n的整数因子。如果n可以分解成整数积f1f2…fm，然后在O(nfi)运算中计算DFT。对于基数为2 的FFT，运算次数为O(nlog2 n)。

基于相同的数学定义，所有FFT函数都提供三种类型的转换:正向、逆向和反向。反向傅里叶变换x = FFT(z)的定义是，

傅里叶逆变换x = IFFT(z)的定义是

 因子1/n使得它是一个真正的逆矩阵。例如，调用gsl_fft_complex_forward()然后调用gsl_fft_complex_inverse()应该返回原始数据(在数值错误内)。

    一般来说，对于变换/逆变换对中的指数符号有两种可能的选择。GSL遵循与FFTPACK相同的约定，使用负指数进行正向转换。这种惯例的优点是，逆变换用简单的傅里叶综合法重建了原来的函数。数值方法使用相反的约定，正指数在正向变换。

反向FFT只是我们对逆FFT的无标度版本的术语，

当结果的整体比例不重要时，使用反向FFT而不是逆通常会很方便，以节省不必要的除法。

16.2 复数FFTs概述

复FFT例程的输入和输出是浮点数的压缩数组。在一个压缩数组中，每个复数的实部和虚部都位于相邻的备用元素中。例如，以下对长度为6的打包数组的定义：

double x[3*2]; gsl_complex_packed_array data = x;

可以通过以下方式用于保存三个复数z [3]的数组：

data[0] = Re(z[0]) data[1] = Im(z[0]) data[2] = Re(z[1]) data[3] = Im(z[1]) data[4] = Re(z[2]) data[5] = Im(z[2])

    数据的数组索引具有与DFT定义相同的顺序，即没有数据的索引转换或排列。

    stride参数允许用户对元素z[stride*i]代替z[i]执行转换。一个大于1的步长可以用来代替矩阵列的位置FFT。步长为1时访问数组，元素之间不需要任何额外的间距。

    要对向量参数（例如gsl_vector_complex * v）执行FFT，请在调用本章中介绍的函数时使用以下定义（或其等效项）：

gsl_complex_packed_array data = v->data; size_t stride = v->stride; size_t n = v->size;

    对于物理应用程序，重要的是要记住，DFT中出现的索引并不直接对应于物理频率。如果DFT的时间步长为∆，则频率域包括正频率和负频率，范围为- 1/(2∆)到0 ~ +1/(2∆)。正频率从数组的开头一直存储到中间，负频率从数组的末尾向后存储。

下面的表格展示了数组数据的布局，以及时域数据z和频域数据x之间的对应关系:

index        z                         x = FFT(z) 0               z(t = 0)               x(f = 0) 1               z(t = 1)              x(f = 1/(n Delta)) 2               z(t = 2)              x(f = 2/(n Delta)) . .               .......                   .................. n/2           z(t = n/2)          x(f = +1/(2 Delta),                                         -1/(2 Delta)) .                ........                  .................. n-3          z(t = n-3)          x(f = -3/(n Delta)) n-2          z(t = n-2)          x(f = -2/(n Delta)) n-1          z(t = n-1)          x(f = -1/(n Delta)

    当n为偶数时，位置n/2包含的正、负频率(+1/(2∆)，-1/(2∆))最大，且相等。如果n是奇数，那么上面表的一般结构仍然适用，但是n/2没有出现。