F # / « accélérateur v2 » implémentation de l'algorithme DFT probablement incorrect

Question

Je suis en train d'expérimenter des concepts de la radio logicielle. De cette article J'ai essayé de mettre en œuvre un GPU-parallélisme transformée de Fourier discrète.

Je suis assez sûr que je pourrais précalculer 90 degrés du péché (i) cos (i), puis il suffit de retourner et répéter plutôt que ce que je fais dans ce code et que ce serait l'accélérer. Mais jusqu'à présent, je ne pense même pas que je obtenir des réponses correctes. Un que j'espère que les contributions donne un résultat 0 tous-zéros, mais toutes les 0,5 entrées en donne 78.9985886f (j'attends 0 résultat dans ce cas aussi). En fait, je suis juste généralement confus. Je n'ai pas de données d'entrée bien et je ne sais pas quoi faire avec le résultat ou la façon de le vérifier.

Cette question est liée à mon autre post

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

ignorer (System.Console.ReadKey ())

Answer 1

Je partage votre grande surprise que votre réponse n'est pas plus proche de zéro. Je vous suggère d'écrire du code naïf pour effectuer votre DFT en F # et de voir si vous pouvez traquer la source de l'écart.

Voici ce que je pense que vous essayez de le faire:

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

Si tout va bien, vous pouvez utiliser ce code naïf pour obtenir une meilleure intuition de la façon dont le code d'accélérateur doit se comporter, ce qui pourrait vous guider dans vos tests du code d'accélérateur. Gardez à l'esprit qu'une partie de la raison de l'écart peut être simplement la précision des calculs - il y a ~ 52 millions d'éléments dans vos tableaux, donc l'accumulation d'une erreur totale de 79 ne peut pas réellement être trop mauvais. FWIW, je reçois une puissance de ~ 0,05 lors de l'exécution du code ci-dessus unique de précision, mais une puissance de ~ 4E-18 lors de l'utilisation du code équivalent à double chiffres de précision.

Answer 2

Deux suggestions:

• vous assurer que vous êtes confus pas en quelque sorte avec degrés radians
• essayer de le faire sans parallélisme, ou tout simplement avec asyncs pour le parallélisme de F #

(En F #, si vous avez un tableau de flotteurs

let a : float[] = ...

alors vous pouvez « ajouter une étape à tous en parallèle » pour produire un nouveau tableau avec

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

(bien que je pense que cela pourrait être plus lent que faire juste une boucle synchrone).)