您的位置:

快速傅里叶变换

一、快速傅里叶变换fft

快速傅里叶变换,即FFT(Fast Fourier Transform),是一种高效的离散傅里叶变换(DFT)算法。FFT能够在O(nlogn)的时间内计算一个长度为n的DFT,相比普通DFT的时间复杂度O(n^2),FFT具有更高的计算效率。


#include <complex>
#include <valarray>
void fft_recursive(std::valarray< std::complex< double > >& x) {
  const size_t N = x.size();
  if (N <= 1) return;
  // divide
  std::valarray< std::complex< double > >& even = x[std::slice(0, N/2, 2)];
  std::valarray< std::complex< double > >& odd = x[std::slice(1, N/2, 2)];
  // conquer
  fft_recursive(even);
  fft_recursive(odd);
  // combine
  for (size_t k = 0; k < N/2; ++k) {
    std::complex< double > t = std::polar(1.0, -2 * M_PI * k / N) * odd[k];
    x[k    ] = even[k] + t;
    x[k+N/2] = even[k] - t;
  }
}

二、快速傅里叶变换蝶形图

快速傅里叶变换的核心思想是蝶形图(Butterfly Diagram),通过对输入序列进行分治和分组,使得序列中的相邻元素具有共性,从而降低计算复杂度。蝶形图是FFT算法的可视化形式,对于长度为N的序列,它可以表示为N/2个蝴蝶,如下图所示:

三、快速傅里叶变换原理

令x(n)为一组离散信号,它的傅里叶变换为X(w),则N个x(n)的DFT可以表示为:

X(k) = ∑n=0N-1 x(n) · e-i2πkn/N

其中e-i2πkn/N为旋转因子,k=0, 1, 2, ..., N-1。FFT利用了DFT的周期性性质,将其分解为多个低阶的DFT,从而降低计算复杂度。

四、快速傅里叶变换标准

FFT算法的标准形式为递归算法和迭代算法。递归算法是一种自上而下的分治算法,通过递归地将问题划分为更小的子问题,并最终将子问题的解组合成原问题的解。迭代算法是一种自下而上的算法,通过循环结构依次对所有蝴蝶进行计算。下面是FFT的迭代算法示例代码:


#include <cmath>
#include <complex>
void fft(std::vector< std::complex< double > >& a) {
  int n = a.size(), L = 31 - __builtin_clz(n);
  static std::vector< std::complex< double > > rt(2, 1);
  for (static int k = 2; k < n; k *= 2) {
    rt.resize(n);
    auto x = std::polar(1.0, M_PI / k);
    for (int i = k; i < 2*k; i++) rt[i] = rt[i-1] * x;
    for (int i = 0; i < n; i += 2 * k) {
      for (int j = 0; j < k; j++) {
        auto z = rt[j + k] * a[i + j + k];
        a[i + j + k] = a[i + j] - z;
        a[i + j] += z;
      }
    }
  }
}

五、快速傅里叶变换

以下是FFT的实现过程:

1、对输入序列进行长度为2的DFT,得到两个长度为1的DFT;

2、对这两个长度为1的DFT进行长度为4的DFT,得到四个长度为1的DFT;

3、对这四个长度为1的DFT进行长度为8的DFT;

4、不断重复,直到得到长度为N的DFT。

六、快速傅里叶变换例题

以下是一个关于FFT的例题,输入一个长度为N的复数序列,输出其DFT的模:


#include <cmath>
#include <iostream>
#include <valarray>
typedef std::complex< double > Complex;
typedef std::valarray< Complex > CArray;
void fft(CArray &x) {
  const size_t N = x.size();
  if (N <= 1) return;
  CArray even = x[std::slice(0, N/2, 2)];
  CArray odd = x[std::slice(1, N/2, 2)];
  fft(even);
  fft(odd);
  for (size_t k = 0; k < N/2; ++k) {
    Complex t = std::polar(1.0, -2 * M_PI * k / N) * odd[k];
    x[k    ] = even[k] + t;
    x[k+N/2] = even[k] - t;
  }
}
int main() {
  int N;
  std::cin >> N;
  CArray data(N);
  for (int i = 0; i < N; ++i) {
    double x, y;
    std::cin >> x >> y;
    data[i] = Complex(x, y);
  }
  fft(data);
  double ans = 0;
  for (auto x : data) ans += std::abs(x);
  std::cout << ans << std::endl;
  return 0;
}

七、快速傅里叶算法

FFT的算法复杂度为O(nlogn),而最快的DFT算法Cooley-Tukey算法也是O(nlogn)。因此,FFT已经成为当前最快的计算DFT的算法之一,广泛应用于音频处理、数字图像处理、信号处理等领域。

八、快速傅里叶变换答案

FFT的输出结果通常是一组复数,其中每个复数的模表示该频率上的振幅,相位表示该频率上的相位。在对离散信号进行频域分析时,我们通常只关注模的值,而忽略相位。

九、快速傅里叶变换算法出自

FFT算法最初由James Cooley和John Tukey在1965年提出,在随后的十几年中,出现了多种针对不同应用场景的FFT变体,这些变体主要是为了提高FFT的计算效率或适应新的信号处理场景。目前,常用的FFT库包括FFTW、MKL和cuFFT等。

十、快速傅里叶变换条件选取

在使用FFT算法时,需要满足以下条件:

1、长度为N的输入序列必须是2的幂次方;

2、如果输入序列不足2的幂次方,需要填充0;

3、如果需要计算逆FFT,需要进行标准化处理;

4、对于实数序列,利用FFT的对称性可以只计算前一半序列,从而降低计算复杂度。