F#/"Accelerator v2" DFT 算法实现可能不正确答案

【问题标题】：F#/"Accelerator v2" DFT algorithm implementation probably incorrectF#/"Accelerator v2" DFT 算法实现可能不正确
【发布时间】：2011-03-08 21:22:30
【问题描述】：

我正在尝试使用软件定义的无线电概念。从article 开始，我尝试实现 GPU 并行离散傅里叶变换。

我很确定我可以预先计算 sin(i) cos(i) 的 90 度，然后只是翻转并重复，而不是我在这段代码中所做的，那样会加快速度。但到目前为止，我什至认为我没有得到正确的答案。如我所料，全零输入给出 0 结果，但作为输入的所有 0.5 给出 78.9985886f （在这种情况下我也预计结果为 0）。基本上，我只是普遍感到困惑。我没有任何好的输入数据，我不知道如何处理结果或如何验证它。

这个问题和我的另一篇帖子here有关

open Microsoft.ParallelArrays
open System

 // X64MulticoreTarget is faster on my machine, unexpectedly
let target = new DX9Target() // new X64MulticoreTarget()

ignore(target.ToArray1D(new FloatParallelArray([| 0.0f |]))) // Dummy operation to warm up the GPU

let stopwatch = new System.Diagnostics.Stopwatch() // For benchmarking

let Hz = 50.0f
let fStep = (2.0f * float32(Math.PI)) / Hz
let shift = 0.0f // offset, once we have to adjust for the last batch of samples of a stream

// If I knew that the periodic function is periodic 
// at whole-number intervals, I think I could keep 
// shift within a smaller range to support streams 
// without overflowing shift - but I haven't 
// figured that out

//let elements = 8192 // maximum for a 1D array - makes sense as 2^13
//let elements = 7240 // maximum on my machine for a 2D array, but why?
let elements = 7240

// need good data!!
let buffer : float32[,] = Array2D.init<float32> elements elements (fun i j -> 0.5f) //(float32(i * elements) + float32(j))) 

let input = new FloatParallelArray(buffer)
let seqN : float32[,] = Array2D.init<float32> elements elements (fun i j -> (float32(i * elements) + float32(j)))
let steps = new FloatParallelArray(seqN)
let shiftedSteps = ParallelArrays.Add(shift, steps)
let increments = ParallelArrays.Multiply(fStep, steps)
let cos_i = ParallelArrays.Cos(increments) // Real component series
let sin_i = ParallelArrays.Sin(increments) // Imaginary component series

stopwatch.Start()
// From the documentation, I think ParallelArrays.Multiply does standard element by 
// element multiplication, not matrix multiplication
// Then we sum each element for each complex component (I don't understand the relationship 
// of this, or the importance of the generalization to complex numbers)
let real = target.ToArray1D(ParallelArrays.Sum(ParallelArrays.Multiply(input, cos_i))).[0]
let imag = target.ToArray1D(ParallelArrays.Sum(ParallelArrays.Multiply(input, sin_i))).[0]
printf "%A in " ((real * real) + (imag * imag)) // sum the squares for the presence of the frequency
stopwatch.Stop()

printfn "%A" stopwatch.ElapsedMilliseconds

忽略 (System.Console.ReadKey())

【问题讨论】：

你得到正确的答案没有并行性？
我不知道如何让它以这种方式运行——我认为它需要一个不同的算法。我仍然不知道正确的答案。

标签： f# radio gpu fft accelerator

【解决方案1】：

我和你一样惊讶，你的答案并不接近于零。我建议编写简单的代码以在 F# 中执行 DFT，并查看是否可以追踪差异的来源。

这就是我认为您正在尝试做的事情：

let N = 7240
let F = 1.0f/50.0f
let pi = single System.Math.PI

let signal = [| for i in 1 .. N*N -> 0.5f |]

let real = 
  seq { for i in 0 .. N*N-1 -> signal.[i] * (cos (2.0f * pi * F * (single i))) }
  |> Seq.sum

let img = 
  seq { for i in 0 .. N*N-1 -> signal.[i] * (sin (2.0f * pi * F * (single i))) }
  |> Seq.sum

let power = real*real + img*img

希望您可以使用这个幼稚的代码来更好地了解加速器代码的行为方式，这可以指导您测试加速器代码。请记住，造成差异的部分原因可能仅仅是计算的精度——数组中有大约 5200 万个元素，因此累积 79 的总误差实际上可能并不算太糟糕。 FWIW，在运行上述单精度代码时，我得到了 ~0.05 的幂，但在使用具有双精度数字的等效代码时，我得到了 ~4e-18 的幂。

【讨论】：

【解决方案2】：

两个建议：

确保您不会以某种方式将度数与弧度混淆
尝试无并行，或仅使用 F# 的异步实现并行

（在 F# 中，如果你有一个浮点数组

let a : float[] = ...

然后您可以“为所有这些并行添加一个步骤”以生成一个新数组

let aShift = a |> (fun x -> async { return x + shift }) 
               |> Async.Parallel |> Async.RunSynchronously

（虽然我预计这可能比只执行同步循环要慢）。

【讨论】：