Fast Fourier Transform

[Buey]

If there's anyone else doing signal processing or the like in AutoIt, here you go. This is for the folks that search the examples forums for existing scripts before hunkering down and writing it themselves. It's even got the syntax for progress bar updating. If it's useful, someone put it in a UDF.

FFT works on input sizes that are powers of 2, so if your input array is not a power of 2, the array will be zero padded to the next power of 2.

Source array format is the same as output array:

Array Format:

x[0] Real Values

x[1] Imaginary Values

x[2] Magnitude

x[3] Phase

Magnitude and phase can be edited out of function, as they are not part of FFT calculation. They are useful values for analysis, which is why I included them.

Usage:

Dim $inputarray[4][500]
LoadSomeValuesIntoArray($inputarray); <-- load some input values into the input array, however you want
$outputarray = _FFT($inputarray); <-- compute FFT over input array.  Output array dimensions: [4][512] (512 is next power of 2 over 500)
$inputarray2 = _FFT($outputarray,-1); <-- compute inverse FFT over output array.  Should arrive with same data that is in inputarray

expandcollapse popup

Func _FFT($srcarray, $dir = 1)
;Fast Fourier Transform: dir 1 for transform, dir -1 for inverse transform
;FFT acts on inputs that are a power of 2 so the function  zero pads the input
;Array Format:
;x[0] Real Values
;x[1] Imaginary Values
;x[2] Magnitude
;x[3] Phase 
    Local $i, $i1, $j, $k, $i2, $l, $l1, $l2, $c1, $c2, $tx, $ty, $t1, $t2, $u1, $u2, $z
    Local $x = $srcarray
    $n = UBound($x,2)   
    $m = Ceiling(Log($n)/Log(2))
    If BitAND($n,$n-1) <> 0 Then; check if power of 2       
        $l = $n
        $n = 2^$m
        ReDim $x[4][$n]
        For $i=$l To $n-1 Step 1
            $x[0][$i] = 0
            $x[1][$i] = 0
            $x[2][$i] = 0
            $x[3][$i] = 0
        Next                
    EndIf
    $i2 = BitShift($n, 1)
    $j = 0
    GUICtrlSetData($progress,0)
    GUICtrlSetData($progresslabel,"Progress: 0%")
    $count = 0
    For $i=0 To $n-2 Step 1
        If $i < $j Then
            $tx = $x[0][$i]
            $ty = $x[1][$i]
            $x[0][$i] = $x[0][$j]
            $x[1][$i] = $x[1][$j]
            $x[0][$j] = $tx
            $x[1][$j] = $ty
        EndIf
        $k = $i2
        While $k <= $j
            $j -= $k
            $k = BitShift($k,1)
        WEnd
        $j += $k
    Next
    $c1 = -1
    $c2 = 0
    $l2 = 1
    For $l=0 To $m-1 Step 1
        $l1 = $l2
        $l2 = BitShift($l2,-1)
        $u1 = 1
        $u2 = 0
        For $j=0 To $l1-1 Step 1
            For $i=$j To $n-1 Step $l2
                $i1 = $i + $l1
                $t1 = $u1 * $x[0][$i1] - $u2 * $x[1][$i1]
                $t2 = $u1 * $x[1][$i1] + $u2 * $x[0][$i1]
                $x[0][$i1] = $x[0][$i] - $t1
                $x[1][$i1] = $x[1][$i] - $t2
                $x[0][$i] += $t1
                $x[1][$i] += $t2
                $count += 1
                GUICtrlSetData($progress,Round(($count)/($m*$n/2)*100))
                GUICtrlSetData($progresslabel,"Progress: " & Round(($count)/($m*$n/2)*100) & "%")
            Next
            $z = $u1 * $c1 - $u2 * $c2
            $u2 = $u1 * $c2 + $u2 * $c1
            $u1 = $z
        Next
        $c2 = Sqrt((1-$c1)/2)
        If $dir > 0 Then
            $c2 *= -1
        EndIf
        $c1 = Sqrt((1+$c1)/2)
    Next
    If $dir < 0 Then
        For $i=0 To $n-1 Step 1
            $x[0][$i] /= $n
            $x[1][$i] /= $n
        Next
    EndIf
    For $i=0 To $n-1 Step 1
        $x[2][$i] = Sqrt($x[0][$i]^2+$x[1][$i]^2)
        $x[3][$i] = ATan($x[1][$i]/$x[0][$i])
    Next
    Return $x
EndFunc

Edited September 13, 2007 by Buey

[Toady]

Cool script. Can you post an example script. What does the scrarray array consist of?

[Buey]

added to original post

[v7ilnytskyy]

On 9/13/2007 at 8:26 PM, Buey said:

Output array Format:

x[0] Real Values

x[1] Imaginary Values

x[2] Magnitude

x[3] Phase

Phase has values: -1, 0, 1.

Can I see the spectrum (frequencies of harmonics and their amplitudes) or somehow calculate their values?