功能惯用的 FFT

Question

我编写这个 radix-2 FFT 的目的是在不牺牲太多性能的情况下使其在功能上符合习惯：

let reverse x bits =
    let rec reverse' x bits y =
        match bits with
        | 0 -> y
        | _ -> ((y <<< 1) ||| (x &&& 1))
               |> reverse' (x >>> 1) (bits - 1) 
     reverse' x bits 0 

let radix2 (vector: Complex[]) (direction: int) =
    let z = vector.Length
    let depth = floor(Math.Log(double z, 2.0)) |> int
    if (1 <<< depth) <> z then failwith "Vector length is not a power of 2"

    // Complex roots of unity; "twiddle factors"
    let unity: Complex[] =
        let xpn = float direction * Math.PI / double z
        Array.Parallel.init<Complex> (z/2) (fun i ->
            Complex.FromPolarCoordinates(1.0, (float i) * xpn))

    // Permutes elements of input vector via bit-reversal permutation
    let pvec = Array.Parallel.init z (fun i -> vector.[reverse i depth])

    let outerLoop (vec: Complex[]) =
        let rec recLoop size =
            if size <= z then
                let mid, step = size / 2, z / size
                let rec inrecLoop i =
                    if i < z then
                        let rec bottomLoop idx k =
                            if idx < i + mid then
                                let temp = vec.[idx + mid] * unity.[k]
                                vec.[idx + mid] <- (vec.[idx] - temp)
                                vec.[idx] <- (vec.[idx] + temp)
                                bottomLoop (idx + 1) (k + step)
                        bottomLoop i 0
                        inrecLoop (i + size)
                inrecLoop 0
                recLoop (size * 2)
        recLoop 2
        vec

    outerLoop pvec

outerLoop 段是我写过的最大的嵌套尾递归混乱。我在 Cooley-Tukey algorithm 的维基百科文章中复制了该算法，但我认为使用高阶函数实现的唯一函数构造会对性能和内存效率造成巨大影响。是否有其他解决方案可以产生相同的结果而不会导致大幅减速，同时仍然是惯用的？

Answer 1

我不是算法工作原理方面的专家，因此可能有一个很好的功能实现，但值得注意的是，使用局部突变在 F# 中是完全惯用的。

你的 radix2 函数从外部是起作用的 - 它需要 vector 数组作为输入，从不改变它，创建一个新数组 pvec 然后初始化（使用一些沿途突变）然后 returns 它。这与 Array.map 等内置函数使用的模式类似（初始化一个新数组，对其进行变异，然后 returns 它）。这通常是一种明智的做事方式，因为一些算法最好使用变异来编写。

在这种情况下，也使用局部可变变量和循环是完全合理的。与尾递归版本相比，这样做将使您的代码更具可读性。我没有测试过这个，但我对你的 outerLoop 函数的天真翻译只是使用三个嵌套循环 - 像这样：

let mutable size = 2
while size <= z do
    let mid, step = size / 2, z / size
    let mutable i = 0
    while i < z do
        for j in 0 .. mid - 1 do
            let idx, k = i + j, step * j
            let temp = pvec.[idx + mid] * unity.[k]
            pvec.[idx + mid] <- (pvec.[idx] - temp)
            pvec.[idx] <- (pvec.[idx] + temp)
        i <- i + size
    size <- size * 2

这可能不完全正确（我这样做只是为了重构您的代码），但我认为它实际上比在这种情况下使用复杂的嵌套尾递归函数更符合习惯。

功能惯用的 FFT

Functionally idiomatic FFT

f#

functional-programming

fft