🌊 快速傅里叶变换 (FFT)

class FFT {
    constructor() {
        this.sinTable = [];
        this.cosTable = [];
    }

    fft(x) {
        const n = x.length;
        if (n === 1) return x;
        if (n % 2 !== 0) throw new Error("序列长度必须是 2 的幂");

        const even = this.fft(x.filter((_, i) => i % 2 === 0));
        const odd = this.fft(x.filter((_, i) => i % 2 === 1));
        const X = new Array(n).fill(0);

        for (let k = 0; k < n / 2; k++) {
            const angle = -2 * Math.PI * k / n;
            const w = {
                real: Math.cos(angle),
                imag: Math.sin(angle)
            };
            const oddTerm = {
                real: w.real * odd[k].real - w.imag * odd[k].imag,
                imag: w.real * odd[k].imag + w.imag * odd[k].real
            };
            X[k] = {
                real: even[k].real + oddTerm.real,
                imag: even[k].imag + oddTerm.imag
            };
            X[k + n / 2] = {
                real: even[k].real - oddTerm.real,
                imag: even[k].imag - oddTerm.imag
            };
        }
        return X;
    }

    fftIterative(x) {
        const n = x.length;
        if ((n & (n - 1)) !== 0) throw new Error("序列长度必须是 2 的幂");

        const X = x.map(c => ({ real: c.real !== undefined ? c.real : c, imag: 0 }));

        // 位逆序重排
        let j = 0;
        for (let i = 1; i < n; i++) {
            let bit = n >> 1;
            while (j & bit) {
                j ^= bit;
                bit >>= 1;
            }
            j ^= bit;
            if (i < j) {
                [X[i], X[j]] = [X[j], X[i]];
            }
        }

        // 逐级蝶形运算
        for (let length = 2; length <= n; length *= 2) {
            const half = length / 2;
            for (let i = 0; i < n; i += length) {
                for (let j = 0; j < half; j++) {
                    const angle = -2 * Math.PI * j / length;
                    const w = {
                        real: Math.cos(angle),
                        imag: Math.sin(angle)
                    };
                    const u = X[i + j];
                    const v = {
                        real: w.real * X[i + j + half].real - w.imag * X[i + j + half].imag,
                        imag: w.real * X[i + j + half].imag + w.imag * X[i + j + half].real
                    };
                    X[i + j] = { real: u.real + v.real, imag: u.imag + v.imag };
                    X[i + j + half] = { real: u.real - v.real, imag: u.imag - v.imag };
                }
            }
        }
        return X;
    }

    magnitudeSpectrum(X) {
        return X.map(c => Math.sqrt(c.real * c.real + c.imag * c.imag));
    }
}

// 测试
const fft = new FFT();
const signal = [];
for (let n = 0; n < 256; n++) {
    const t = n / 1000;
    signal.push({
        real: Math.sin(2 * Math.PI * 50 * t) + 0.5 * Math.sin(2 * Math.PI * 120 * t),
        imag: 0
    });
}
const X = fft.fftIterative(signal);
const mags = fft.magnitudeSpectrum(X);
console.log("峰值频率:", mags.indexOf(Math.max(...mags.slice(1, 128))) * 1000 / 256, "Hz");

package main

import (
    "fmt"
    "math"
)

type Complex struct {
    real, imag float64
}

func (c Complex) Add(o Complex) Complex {
    return Complex{c.real + o.real, c.imag + o.imag}
}

func (c Complex) Sub(o Complex) Complex {
    return Complex{c.real - o.real, c.imag - o.imag}
}

func (c Complex) Mul(o Complex) Complex {
    return Complex{
        c.real*o.real - c.imag*o.imag,
        c.real*o.imag + c.imag*o.real,
    }
}

func FFT(x []Complex) []Complex {
    n := len(x)
    if n == 1 {
        return x
    }

    var even, odd []Complex
    for i := 0; i < n; i++ {
        if i%2 == 0 {
            even = append(even, x[i])
        } else {
            odd = append(odd, x[i])
        }
    }

    evenFFT := FFT(even)
    oddFFT := FFT(odd)

    X := make([]Complex, n)
    for k := 0; k < n/2; k++ {
        angle := -2 * math.Pi * float64(k) / float64(n)
        w := Complex{math.Cos(angle), math.Sin(angle)}
        X[k] = evenFFT[k].Add(w.Mul(oddFFT[k]))
        X[k+n/2] = evenFFT[k].Sub(w.Mul(oddFFT[k]))
    }
    return X
}

func main() {
    n := 256
    var signal []Complex
    for i := 0; i < n; i++ {
        t := float64(i) / 1000.0
        val := math.Sin(2*math.Pi*50*t) + 0.5*math.Sin(2*math.Pi*120*t)
        signal = append(signal, Complex{val, 0})
    }

    X := FFT(signal)
    
    var maxMag float64
    var peakIdx int
    for i := 1; i < n/2; i++ {
        mag := math.Sqrt(X[i].real*X[i].real + X[i].imag*X[i].imag)
        if mag > maxMag {
            maxMag = mag
            peakIdx = i
        }
    }

    freqRes := 1000.0 / float64(n)
    fmt.Printf("频率分辨率：%.2f Hz\n", freqRes)
    fmt.Printf("峰值频率：%.1f Hz\n", float64(peakIdx)*freqRes)
}