PureBoard

Code: Alles auswählen

EnableExplicit

Structure Complex
	r.d
	i.d
EndStructure
Procedure DFT_2D(*input.Double, *output.Complex, width.i, height.i)
	Protected i.i, j.i, x.i, y.i, d.d, *pixel.Double, *pixelOut.Complex
	Protected invSqSize.d = 1. / Sqr(width * height)
	*pixelOut = *output
	For i = 0 To width - 1
		For j = 0 To height - 1
			*pixel = *input
			*pixelOut\r = 0.0
			*pixelOut\i = 0.0
			For x = 0 To width - 1
				For y = 0 To height -1
					;exp(-imag * 2 * pi * i * x / width) * exp(-imag * 2 * pi * j * y / height)
					;= cos((2 * pi * x * i) / width + (2 * pi * y * j) / height) - imag * sin((2 * pi * x * i) / width + (2 * pi * y * j) / height)
					d = 2 * #PI * ((x * i) / width + (y * j) / height)
					*pixelOut\r + (*pixel\d * Cos(d))
					*pixelOut\i + (*pixel\d * -Sin(d))
					*pixel + SizeOf(Double)
				Next
			Next
			;Debug StrD(*pixelOut\r, 4) + " + " + StrD(*pixelOut\i, 4) + "i"
			*pixelOut\r * invSqSize
			*pixelOut\i * invSqSize
			*pixelOut + SizeOf(Complex)
		Next
	Next
EndProcedure

Procedure InvDFT_2D(*input.Complex, *output.Double, width.i, height.i)
	Protected i.i, j.i, x.i, y.i, *pixel.Complex, *pixelOut.Double
	Protected real.d, imag.d, d.d, tmpReal.d, tmpImag.d
	Protected invSqSize.d = 1. / Sqr(width * height)
	*pixelOut = *output
	For i = 0 To width - 1
		For j = 0 To height - 1
			*pixel = *input
			tmpReal = 0.0
			tmpImag = 0.0
			For x = 0 To width - 1
				For y = 0 To height -1
					;exp(imag * 2 * pi * i * x / width) * exp(imag * 2 * pi * j * y / height)
					;= cos((2 * pi * x * i) / width -(2 * pi * y * j) / height) - imag * sin((2 * pi * x * i) / width - (2 * pi * y * j) / height)
					d = (2 * #PI * x * i) / width + (2 * #PI * y * j) / height
					real = Cos(d)
					imag = Sin(d)
					tmpReal + (real * *pixel\r - imag * *pixel\i)
					tmpImag + (real * *pixel\i + imag * *pixel\r)
					*pixel + SizeOf(Complex)
				Next
			Next
			*pixelOut\d = tmpReal * invSqSize
			*pixelOut + SizeOf(Double)
		Next
	Next
EndProcedure

#width = 2
#height = 2

Dim image.d(#width - 1, #height - 1)
Dim fft.Complex(#width - 1, #height - 1)
Dim image2.d(#width - 1, #height - 1)

Define timeGenerate.i, timeDFT.i, timeInvDFT.i
Define i.i, j.i



timeGenerate = ElapsedMilliseconds()
For i = 0 To #width - 1
	For j = 0 To #height - 1
		image(i, j) = Random(100) / 100.
	Next
Next
timeGenerate = ElapsedMilliseconds() - timeGenerate
; image(0, 0) = 4			;7.07
; image(0, 1) = 3.5		;1.06 - 2.06i
; image(0, 2) = 4			;0
; image(0, 3) = 3.5		;1.06 + 0.06i
; image(0, 4) = 1			;0
; image(0, 5) = 1.5		;1.06 - 0.06i
; image(0, 6) = 1			;0
; image(0, 7) = 1.5		;1.06 + 2.06i

Debug "-------"
timeDFT = ElapsedMilliseconds()
DFT_2D(@image(), @fft(), #width, #height)
timeDFT = ElapsedMilliseconds() - timeDFT

For i = 0 To #width - 1
	For j = 0 To #height - 1
		Debug StrD(image(i, j), 4)
	Next
Next

Debug "-------"

For i = 0 To #width - 1
	For j = 0 To #height - 1
		Debug StrD(fft(i, j)\r, 4) + " + " + StrD(fft(i, j)\i, 4) + "i"
	Next
Next

Debug "-------"
timeInvDFT = ElapsedMilliseconds()
InvDFT_2D(@fft(), @image2(), #width, #height)
timeInvDFT = ElapsedMilliseconds() - timeInvDFT

For i = 0 To #width - 1
	For j = 0 To #height - 1
		Debug StrD(image2(i, j), 4)
	Next
Next

MessageRequester("Running time", "Generate: " + Str(timeGenerate) + " ms" + #LF$ + "DFT: " + Str(timeDFT) + " ms" + #LF$ + "Inverse DFT: " + Str(timeInvDFT) + " ms")