傅里叶变换：频域的世界

时域的卷积 = 频域的乘法。这个定理是信号处理最重要的计算加速工具。

离散傅里叶变换（DFT）

DFT 把 N 个时域采样变换为 N 个频域系数：X[k] = Σ_{n=0}^{N-1} x[n] × e^{-2πi·kn/N}。朴素实现是 O(N²)，因为每个频率系数需要与所有 N 个时域样本做点积。

package main

import (
    "fmt"
    "math"
    "math/cmplx"
)

// 朴素 DFT：O(N²)
func dft(x []float64) []complex128 {
    N := len(x)
    X := make([]complex128, N)
    for k := 0; k < N; k++ {
        for n := 0; n < N; n++ {
            angle := -2 * math.Pi * float64(k*n) / float64(N)
            X[k] += complex(x[n], 0) * cmplx.Exp(complex(0, angle))
        }
    }
    return X
}

// Cooley-Tukey FFT：O(N log N)，N 必须是 2 的幂
func fft(x []complex128) []complex128 {
    N := len(x)
    if N <= 1 {
        return x
    }
    // 分治：偶数下标和奇数下标
    even := make([]complex128, N/2)
    odd := make([]complex128, N/2)
    for i := 0; i < N/2; i++ {
        even[i] = x[2*i]
        odd[i] = x[2*i+1]
    }
    E := fft(even)
    O := fft(odd)
    result := make([]complex128, N)
    for k := 0; k < N/2; k++ {
        t := cmplx.Exp(complex(0, -2*math.Pi*float64(k)/float64(N))) * O[k]
        result[k] = E[k] + t
        result[k+N/2] = E[k] - t
    }
    return result
}

func main() {
    // 合成信号：5Hz + 10Hz 正弦波
    N := 64
    signal := make([]float64, N)
    for i := range signal {
        t := float64(i) / float64(N)
        signal[i] = math.Sin(2*math.Pi*5*t) + 0.5*math.Sin(2*math.Pi*10*t)
    }

    // FFT
    cx := make([]complex128, N)
    for i, v := range signal { cx[i] = complex(v, 0) }
    X := fft(cx)

    // 找主要频率分量
    fmt.Println("主要频率分量（振幅 > 0.1N）:")
    for k := 0; k < N/2; k++ {
        amp := cmplx.Abs(X[k])
        if amp > 0.1*float64(N) {
            fmt.Printf("  频率 %2d Hz: 振幅 = %.1f\n", k, amp/float64(N)*2)
        }
    }
}

卷积定理与音频处理

时域卷积 (f * g)[n] = Σ f[m]g[n-m] 的复杂度是 O(N²)，但在频域只是点乘：FFT(f * g) = FFT(f) × FFT(g)。所以大卷积的算法是：FFT → 点乘 → IFFT，复杂度降到 O(N log N)。音频均衡器、图像模糊、CNN 的卷积层都用这个技巧。

package main

import (
    "fmt"
    "math"
    "math/cmplx"
)

// 用 FFT 做快速多项式乘法（卷积定理的应用）
// a = [1,2,3]，b = [4,5,6]，求 a*b 展开系数
func polyMul(a, b []float64) []float64 {
    // 补零到 2 的幂次
    n := 1
    for n < len(a)+len(b) { n <<= 1 }

    ca := make([]complex128, n)
    cb := make([]complex128, n)
    for i, v := range a { ca[i] = complex(v, 0) }
    for i, v := range b { cb[i] = complex(v, 0) }

    fa := fft(ca)
    fb := fft(cb)

    // 频域点乘
    fc := make([]complex128, n)
    for i := range fc { fc[i] = fa[i] * fb[i] }

    // IFFT（反转后 FFT 再除以 n）
    for i := 1; i < n/2; i++ { fc[i], fc[n-i] = fc[n-i], fc[i] }
    result := fft(fc)
    out := make([]float64, len(a)+len(b)-1)
    for i := range out {
        out[i] = math.Round(real(result[i])/float64(n)*10) / 10
    }
    return out
}

func fft(x []complex128) []complex128 {
    N := len(x)
    if N <= 1 { return x }
    even, odd := make([]complex128, N/2), make([]complex128, N/2)
    for i := 0; i < N/2; i++ { even[i] = x[2*i]; odd[i] = x[2*i+1] }
    E, O := fft(even), fft(odd)
    res := make([]complex128, N)
    for k := 0; k < N/2; k++ {
        t := cmplx.Exp(complex(0, -2*math.Pi*float64(k)/float64(N))) * O[k]
        res[k] = E[k] + t
        res[k+N/2] = E[k] - t
    }
    return res
}

func main() {
    // (1 + 2x + 3x²)(4 + 5x + 6x²) = 4 + 13x + 28x² + 27x³ + 18x⁴
    result := polyMul([]float64{1, 2, 3}, []float64{4, 5, 6})
    fmt.Println(result)  // [4 13 28 27 18]
}

傅里叶的洞见：任何信号都是正弦波的叠加。这一句话改变了整个现代通信。

— The Fourier Transform and Its Applications