A frequency is a question you ask the signal频率,是向信号提出的一个问题

Suppose eight sensors record a repeating motion. The DFT asks how strongly those eight numbers match each of eight rotating patterns. Every output is a weighted sum of the same input samples. FFT is a way to reuse those sums; it computes the same transform, not an approximation with fewer frequencies.假设八个采样点记录了一段周期运动。DFT 会询问:这八个数分别与八种旋转模式有多匹配?每个输出都是对同一组输入做加权求和。FFT 通过复用这些求和来加速,计算的仍然是同一个变换,不是少看几个频率的近似。

1.1

Here i is the imaginary unit, the forward transform uses a negative exponent, and there is no 1/n factor. The inverse uses the opposite sign and divides by n. Fix these conventions before comparing implementations: a sign or normalization mismatch can look like an algorithm bug.这里 i 是虚数单位,正向变换采用负指数,不带 1/n 因子。逆变换采用相反符号并除以 n。比较两个实现前,要先对齐这些约定,否则符号或归一化差异会被误判成算法错误。

Change the signal, watch the spectrum换一种信号,看频谱如何变化

Choose an eight-sample signal. The table shows the signed samples; the bars show unnormalized magnitudes |Xₖ|. Frequency k represents k/8 cycles per sample. The browser computes this tiny display with the direct DFT so it can serve as a transparent reference for the FFT code below.选择一组八点信号。表格显示带正负号的采样值,柱状图显示未归一化的幅值 |Xₖ|。第 k 个频率对应每个采样间隔 k/8 个周期。浏览器用直接 DFT 计算这个小图,便于对照下方 FFT 代码。

Magnitude spectrum · shared numerical results幅度谱 · 中英文共用同一组结果

An impulse at j = 0 contributes 1 to every frequency. A constant signal adds coherently only at k = 0, giving 8. Alternating +1 and −1 concentrates at k = 4. A unit sine at bin 1 gives magnitudes 4 at bins 1 and 7: the second peak is its negative-frequency partner, not a second unrelated oscillation.j = 0 处的单位脉冲对每个频率都贡献 1。常量信号只在 k = 0 同相累加,得到 8。正负交替的信号集中在 k = 4。索引为 1 的单位正弦,在 1 和 7 处各有幅值 4:后者是负频率伙伴,不是另一个无关振动。

Split even and odd positions把偶数位置和奇数位置分开

Let n = 2m. Group the terms whose input index is 2r and those whose index is 2r+1. The even-index sum becomes an m-point DFT, as does the odd-index sum. Denote them E and O. One rotation of O can be reused to produce two outputs.令 n = 2m,把输入索引为 2r 和 2r+1 的项分组。偶数位置的求和变成 m 点 DFT,奇数位置也一样,分别记为 E 和 O。把 O 旋转一次,就能同时得到两个输出。

3.1
3.2

Why does the lower output use a minus sign? Because ωₙᵐ = −1, so increasing the output index by m flips the odd contribution and leaves the even contribution unchanged. This two-input, two-output pattern is the butterfly.为什么后一半要用减号?因为 ωₙᵐ = −1,输出索引增加 m 后,奇数位置的贡献翻转符号,偶数位置的贡献不变。这个两路输入、两路输出的组合就是蝶形运算。

Eₖ+ tXₖ
Oₖ → t = ωₙᵏOₖ Eₖ − tXₖ₊ₘ

For x = [1,2,3,4], the even part [1,3] transforms to [4,−2] and the odd part [2,4] to [6,−2]. With twiddles [1,−i], merging gives [10,−2+2i,−2,−2−2i]. Check this four-point case by hand before thinking about thousands of samples.取 x = [1,2,3,4],偶数位置 [1,3] 的变换为 [4,−2],奇数位置 [2,4] 为 [6,−2]。乘上旋转因子 [1,−i] 再合并,得到 [10,−2+2i,−2,−2−2i]。先手算通这个四点例子,再考虑成千上万个采样点。

Where the speedup comes from加速来自哪里?

4.1

At every level, the total merge work is linear in n, and there are log₂n levels. For n = 1024, direct evaluation has 1,048,576 input–frequency terms; radix-2 merging has 5,120 butterflies. These are different kinds of work, so their ratio is not a measured speedup or an exact FLOP ratio.每一层的合并总工作量与 n 成正比,共有 log₂n 层。n = 1024 时,直接求和有 1,048,576 个输入与频率的配对项;基 2 合并共有 5,120 个蝶形。两种计数单位不同,因此它们的比值既不是实测加速比,也不是精确的 FLOP 比值。

This recursive version splits and allocates temporary vectors to make the idea visible. It uses O(n) peak auxiliary storage, but allocates and copies across recursion levels. Iterative implementations can work in place after a suitable permutation; tuned libraries also plan transform decompositions and memory access.这个递归版本通过拆分与临时向量把结构写清楚,峰值辅助空间为 O(n),但会跨递归层分配和复制。迭代实现可以在适当重排后原地计算;经过调优的库还会规划变换分解方式和访存过程。

Run FFT against a direct DFT用直接 DFT 校验 FFT

#include <algorithm>
#include <cmath>
#include <complex>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <vector>
using C = std::complex<double>;
using Signal = std::vector<C>;
const double pi=std::acos(-1.0);

// Forward transform: negative exponent, no normalization.
Signal radix2(const Signal& x) {
    const std::size_t n=x.size();
    if (n==0 || (n&(n-1))!=0)
        throw std::invalid_argument("size must be a positive power of two");
    if (n==1) return x;
    Signal even(n/2),odd(n/2),y(n);
    for (std::size_t j=0;j<n/2;++j) {even[j]=x[2*j];odd[j]=x[2*j+1];}
    auto e=radix2(even),o=radix2(odd);
    for (std::size_t k=0;k<n/2;++k) {
        const C t=std::polar(1.0,-2*pi*double(k)/double(n))*o[k];
        y[k]=e[k]+t;
        y[k+n/2]=e[k]-t;
    }
    return y;
}
Signal direct_dft(const Signal& x) {
    Signal y(x.size());
    for (std::size_t k=0;k<x.size();++k)
        for (std::size_t j=0;j<x.size();++j)
            y[k]+=x[j]*std::polar(1.0,-2*pi*double(j)*double(k)/double(x.size()));
    return y;
}
Signal inverse_fft(Signal y) {
    for (auto& v:y) v=std::conj(v);
    y=radix2(y);
    for (auto& v:y) v=std::conj(v)/double(y.size());
    return y;
}
double error(const Signal& a,const Signal& b) {
    double result=0;
    for (std::size_t i=0;i<a.size();++i) result=std::max(result,std::abs(a[i]-b[i]));
    return result;
}
int main() {
    Signal x(8);
    for (std::size_t j=0;j<x.size();++j) x[j]=std::sin(2*pi*double(j)/8);
    auto y=radix2(x);
    std::cout << std::fixed << std::setprecision(3) << "magnitudes:";
    for (auto v:y) std::cout << ' ' << std::abs(v);
    std::cout << '\n';
    double dft_error=0,roundtrip_error=0;
    for (int n : {1,2,4,8,16,32}) {
        Signal input(n);
        for (int j=0;j<n;++j) input[j]=C(std::sin(0.7*j),std::cos(0.3*j));
        auto result=radix2(input);
        dft_error=std::max(dft_error,error(result,direct_dft(input)));
        roundtrip_error=std::max(roundtrip_error,error(input,inverse_fft(result)));
    }
    bool rejected=false;
    try {radix2(Signal(3));} catch(const std::invalid_argument&) {rejected=true;}
    std::cout << std::scientific << "DFT error=" << dft_error
              << " roundtrip error=" << roundtrip_error
              << " rejected size 3=" << rejected << '\n';
    return dft_error<1e-10 && roundtrip_error<1e-10 && rejected ? 0 : 1;
}

The program compares both methods for six power-of-two sizes, checks an inverse round trip, and rejects size 3. The sine example should have magnitudes [0,4,0,0,0,0,0,4], up to rounding. Small nonzero errors are expected with floating-point complex arithmetic.程序对六种 2 的幂次长度比较两种算法,检查正逆变换往返,并拒绝长度 3。正弦示例的幅值应接近 [0,4,0,0,0,0,0,4],允许舍入误差。浮点复数运算产生很小的非零误差是正常现象。

Read the picture without overclaiming读懂图,也要读懂边界

Radix-2 requires a positive power-of-two length; FFT as a family does not. Zero-padding a length-6 signal to length 8 changes the sampled frequency grid: it is not the original six-point DFT. A frequency between bins spreads energy across bins because a finite observation window generally cuts an incomplete number of periods.基 2 算法要求长度是正的 2 的幂,但 FFT 这一算法家族并不限于此。把六点信号补零到八点,会改变频率采样网格,不等于原来的六点 DFT。若频率落在两个网格点之间,有限观测窗通常截不到整数个周期,能量就会散布到多个频率点。

Magnitude also discards phase. Two spectra with the same bar heights can reconstruct different time signals. Try changing the C++ input to a shifted impulse: magnitudes remain one, while the complex phases rotate.幅值还丢掉了相位信息。柱子一样高的两组频谱,也可能重建出不同的时域信号。可以把 C++ 输入改为平移后的单位脉冲:幅值仍全为一,复数相位却会旋转。

Read further继续阅读

FFTW: transform sign, normalization, and frequency orderingFFTW:变换符号、归一化和频率排列

Explore phase with a shared wave experiment用波动实验继续理解相位