功能上惯用的FFT

Question

我编写这个基-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

Answer 1

我不是算法工作的专家，所以可能有一个很好的函数实现，但值得注意的是，在F#中使用局部变异是非常惯用的。

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

在这种情况下，使用局部可变变量和循环是完全合理的。与尾递归版本相比，这样做将使您的代码更具可读性。我没有对此进行测试，但我对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

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