F #/"Accelerator v2" Implémentation de l'algorithme DFT probablement incorrecte

Question

J'essaie d'expérimenter des concepts radio définis par logiciel. De ce article j'ai essayé de mettre en application une transformation de Fourier Discrète de GPU-parallélisme. Je suis assez sûr que je pourrais pré-calculer 90 degrés du sin (i) cos (i) et puis juste retourner et répéter plutôt que ce que je fais dans ce code et que cela l'accélérerait. Mais jusqu'à présent, je ne pense même pas que j'obtiens des réponses correctes. Une entrée tout-zéros donne un résultat de 0 comme je m'y attendais, mais tous les 0.5 comme entrées donnent 78.9985886f (je m'attendrais à un résultat de 0 dans ce cas aussi). Fondamentalement, je suis généralement confus. Je n'ai pas de bonnes données d'entrée et je ne sais pas quoi faire avec le résultat ou comment le vérifier.

Cette question est liée à mon autre post 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 

ignorer (System.Console.ReadKey())

Answer 1

Je partage votre surprise que votre réponse ne soit pas plus proche de zéro. Je suggère d'écrire du code naïf pour effectuer votre DFT en F # et voir si vous pouvez trouver la source de la divergence.

Voici ce que je pense que vous essayez de 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 

Espérons que 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 de l'accélérateur. Gardez à l'esprit qu'une partie de la raison de l'écart peut simplement être la précision des calculs - il y a environ 52 millions d'éléments dans vos tableaux, donc accumuler une erreur totale de 79 peut ne pas être vraiment mauvais. FWIW, je reçois une puissance de ~ 0.05 lors de l'exécution du code simple précision ci-dessus, mais une puissance de ~ 4e-18 lors de l'utilisation de code équivalent avec des nombres double précision.

Answer 2

Deux suggestions:

• vous n'êtes pas assurer une certaine façon déroutante degrés avec radians
• essayez de le faire sans-parallélisme, ou juste avec asynchrones F # pour le parallélisme

(In 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 

(même si je pense que cela pourrait être plus lent que de faire une boucle synchrone).)