#COMPILE EXE
#DIM ALL

FUNCTION PBMAIN() AS LONG
    LOCAL data(16) AS DOUBLE
    LOCAL real_out(16) AS DOUBLE
    LOCAL imag_out(16) AS DOUBLE
    LOCAL n AS LONG
    LOCAL i AS LONG
    LOCAL j AS LONG
    LOCAL k AS LONG
    LOCAL len AS LONG
    LOCAL angle AS DOUBLE
    LOCAL cos_val AS DOUBLE
    LOCAL sin_val AS DOUBLE
    LOCAL temp_r AS DOUBLE
    LOCAL temp_i AS DOUBLE
    LOCAL pi_val AS DOUBLE
    LOCAL temp_str AS STRING

    pi_val = 3.14159265

    PRINT "FFT (Fast Fourier Transform) Demo"
    PRINT "================================"
    PRINT ""

    ' Simple sine wave input: sin(2*pi*3*t) with 8 samples
    n = 8
    i = 0
    WHILE i < n
        angle = 2.0 * pi_val * 3.0 * i / n
        data(i) = SIN(angle)
        real_out(i) = data(i)
        imag_out(i) = 0.0
        i = i + 1
    WEND

    PRINT "Input samples:"
    temp_str = ""
    i = 0
    WHILE i < n
        temp_str = temp_str + STR$(data(i))
        IF i < n - 1 THEN temp_str = temp_str + " "
        i = i + 1
    WEND
    PRINT temp_str

    ' FFT: decimation-in-time Cooley-Tukey
    len = n
    j = 0
    i = 0
    WHILE i < n - 1
        IF i < j THEN
            temp_r = real_out(j)
            temp_i = imag_out(j)
            real_out(j) = real_out(i)
            imag_out(j) = imag_out(i)
            real_out(i) = temp_r
            imag_out(i) = temp_i
        END IF
        k = n / 2
        WHILE k >= 1 AND k <= j
            j = j - k
            k = k / 2
        WEND
        j = j + k
        i = i + 1
    WEND

    len = 2
    WHILE len <= n
        angle = 2.0 * pi_val / len
        cos_val = COS(angle)
        sin_val = SIN(angle)

        i = 0
        WHILE i < n
            temp_r = 1.0
            temp_i = 0.0
            j = 0
            WHILE j < len / 2
                k = i + j
                ' butterfly
                temp_r = real_out(k + len / 2)
                temp_i = imag_out(k + len / 2)
                real_out(k + len / 2) = real_out(k) - temp_r
                imag_out(k + len / 2) = imag_out(k) - temp_i
                real_out(k) = real_out(k) + temp_r
                imag_out(k) = imag_out(k) + temp_i
                j = j + 1
            WEND
            i = i + len
        WEND
        len = len * 2
    WEND

    PRINT ""
    PRINT "FFT output (magnitude):"
    i = 0
    WHILE i < n
        temp_str = STR$(SQR(real_out(i) * real_out(i) + imag_out(i) * imag_out(i)))
        PRINT "  bin " + STR$(i) + ": " + temp_str
        i = i + 1
    WEND

    ' Test with simple DC signal
    n = 4
    real_out(0) = 1.0
    imag_out(0) = 0.0
    real_out(1) = 1.0
    imag_out(1) = 0.0
    real_out(2) = 1.0
    imag_out(2) = 0.0
    real_out(3) = 1.0
    imag_out(3) = 0.0

    j = 0
    i = 0
    WHILE i < 3
        IF i < j THEN
            temp_r = real_out(j)
            real_out(j) = real_out(i)
            real_out(i) = temp_r
        END IF
        k = 2
        WHILE k >= 1 AND k <= j
            j = j - k
            k = k / 2
        WEND
        j = j + k
        i = i + 1
    WEND

    len = 2
    WHILE len <= n
        angle = 2.0 * pi_val / len
        cos_val = COS(angle)
        sin_val = SIN(angle)
        i = 0
        WHILE i < n
            j = 0
            WHILE j < len / 2
                k = i + j
                temp_r = real_out(k + len / 2)
                real_out(k + len / 2) = real_out(k) - temp_r
                real_out(k) = real_out(k) + temp_r
                j = j + 1
            WEND
            i = i + len
        WEND
        len = len * 2
    WEND

    PRINT ""
    PRINT "DC signal FFT (all 1s): should show peak at bin 0"
    i = 0
    WHILE i < n
        temp_str = STR$(real_out(i))
        PRINT "  bin " + STR$(i) + ": " + temp_str
        i = i + 1
    WEND

    PRINT ""
    PRINT "FFT demo complete."

    FUNCTION = 0
END FUNCTION
