I would suggest that you open a topic under General and maybe I or someone else can help you.
Currently, I've a COVID infection...
I would suggest that you open a topic under General and maybe I or someone else can help you.
Code: Select all
const pi = 4 * atn(1)
'' ---------------
' -~- FFT -~-
' ----------------
''2 ^ 10 = 1024
''2 ^ 18 = 262144
Dim As uByte lutpow = 14
Dim Shared As Integer MaxPowerOfTwo = 20
If lutpow > MaxPowerOfTwo Then lutpow = MaxPowerOfTwo
Dim As Integer size_PADLUT = 1 Shl lutpow
Function NumberOfBitsNeeded(ByVal PowerOfTwo as Integer) as Integer
for I As Integer= 0 to MaxPowerOfTwo
if (PowerOfTwo and (1 shl I)) <> 0 then
return I
End If
Next
Return 0
End Function
Function ReverseBits(ByVal Index As Integer,ByVal NumBits As Integer) As Integer
Dim As Integer Rev
for I As Integer = 0 to NumBits - 1
Rev Shl= 1
Rev Or= Index And 1
Index Shr= 1
Next
Return Rev
End Function
Dim Shared As Integer RevBits(size_PADLUT)
Dim As Integer BitsNeeded = NumberOfBitsNeeded(size_PADLUT)
For I As Integer = 0 To size_PADLUT
RevBits(I) = ReverseBits(I, BitsNeeded)
Next
Function IsPowerOfTwo(ByVal X as Integer) As Integer
Dim As Integer I, Y
Y = 2
for I = 1 To MaxPowerOfTwo
if X = Y then Return TRUE
Y Shl= 1
Next
Return FALSE
End Function
Sub FourierTransform(ByVal AngleNumerator As Single, _
ByVal NumSamples As Integer, _
lpReal As Single Ptr, _
lpImag As Single Ptr, _
lpRealOut As Single Ptr, _
lpImagOut As Single Ptr)
' http://www.koders.com/delphi/fidB6DD10205DAFD71EFF5093D4916237BB2801DD8D.aspx
Dim As Integer I, J, K, N, BlockSize, BlockEnd', NumBits
Dim As single Delta_angle, Delta_ar
Dim As single Alpha, Beta
Dim As single Tr, Ti, Ar, Ai
if not IsPowerOfTwo(NumSamples) or (NumSamples < 2) then
?"Error in procedure Fourier: NumSamples="; NumSamples
?" is not a positive integer power of 2."
Exit Sub
End If
'NumBits = NumberOfBitsNeeded(NumSamples)
for I = 0 to NumSamples
J = RevBits(I)
lpRealOut[J] = lpReal[I]
lpImagOut[J] = lpImag[I]
Next
BlockEnd = 1
BlockSize = 2
while BlockSize <= NumSamples
Delta_angle = AngleNumerator / BlockSize
Alpha = Sin(0.5 * Delta_angle)
Alpha = 2.0 * Alpha * Alpha
Beta = Sin(Delta_angle)
I = 0
while I < NumSamples
Ar = 1.0 ''cos(0)
Ai = 0.0 ''sin(0)
J = I
for N = 0 to BlockEnd - 1
K = J + BlockEnd
Tr = Ar * lpRealOut[K] - Ai * lpImagOut[K]
Ti = Ar * lpImagOut[K] + Ai * lpRealOut[K]
lpRealOut[K] = lpRealOut[J] - Tr
lpImagOut[K] = lpImagOut[J] - Ti
lpRealOut[J] += Tr
lpImagOut[J] += Ti
Delta_ar = Alpha * Ar + Beta * Ai
Ai -= (Alpha * Ai - Beta * Ar)
Ar -= Delta_ar
J += 1
Next
I += BlockSize
Wend
BlockEnd = BlockSize
BlockSize Shl= 1
Wend
End Sub
type signal_and_transform
declare sub resize( as long )
declare sub FFT( as long )
as single sr(any)
as single si(any)
as single dr(any)
as single di(any)
end type
sub signal_and_transform.resize( c as long )
redim sr(c-1)
redim si(c-1)
redim dr(c-1)
redim di(c-1)
end sub
sub signal_and_transform.FFT( c_samps as long )
FourierTransform( _
2 * PI, c_samps, @sr(0), @si(0), @dr(0), @di(0) )
end sub
dim as signal_and_transform signal
var num_samples = 2 ^ 9
signal.resize num_samples
for i as long = 0 to num_samples - 1
signal.sr(i) = (rnd - .5)
signal.si(i) = (rnd - .5)
next
signal.FFT num_samples
var w = 800
var h = 600
screenres w, h
for i as long = 0 to num_samples - 1
pset ( i, h/2*( 1 + signal.dr(i)/sqr(num_samples) ) )
next
sleep