Author Topic: Fast Fourier Transform  (Read 100 times)

guga

  • Moderator
  • Member
  • *****
  • Posts: 824
  • Assembly is a state of art.
    • RosAsm
Fast Fourier Transform
« on: April 04, 2017, 05:17:53 AM »
Hi guys.

I ported a Fast Fourier Transform used in pHash algo to RosAsm. I hope it can be usefull for someone. It seems slow, though :(

Code: [Select]

[phcomplex.reDis 0
 phcomplex.imDis 8]

[Size_Of_phcomplex 16]

[Float_Two_PI: R$ 6.28318530717958623]
[TmpFloat1: R$ 0]
[Float_One: R$ 1.0]

; InputX: Input array of Real data (Float Double)
; nfactor: The total amount of elements on the array
; XFactor: A buffer that will store the resultant values. (In Double float too.)

Proc fft:
    Arguments @InputX, @nfactor, @XFactor
    Local @twiddle_factors, @XtStruct, @kCounter, @MaxCount, @CurTwiddle_factors
    Uses ecx, edx

    mov eax D@nfactor | shr eax 1 | shl eax 4
    C_call 'msvcrt.malloc' eax
    mov D@twiddle_factors eax
    mov D@CurTwiddle_factors eax

    mov eax D@nfactor | shl eax 4 ; mul by Size_Of_phcomplex
    C_call 'msvcrt.malloc' eax
    mov D@XtStruct eax

    mov D@kCounter 0

    mov ebx D@nfactor | shr ebx 1

    .While D@kCounter < ebx

        fild F@kCounter | fmul R$Float_Two_PI | fidiv F@nfactor | fstp R$TmpFloat1
        call polar_to_complex D@CurTwiddle_factors, Float_One, TmpFloat1

        add D@CurTwiddle_factors Size_Of_phcomplex
        inc D@kCounter
    .End_While

    call fft_calc D@nfactor, D@InputX, D@XFactor, D@XtStruct, 1, D@twiddle_factors
    C_call 'msvcrt.free' D@twiddle_factors
    C_call 'msvcrt.free' D@XtStruct
    xor eax eax

EndP
_________________________________
Code: [Select]
Proc polar_to_complex:
    Arguments @Output, @RValue, @Theta
    Uses edi, ebx

    mov eax D@Output
    mov ebx D@Theta
    mov edi D@RValue
    fld R$ebx | fcos | fmul R$edi | fstp R$eax+phcomplex.reDis
    fld R$ebx | fsin | fmul R$edi | fstp R$eax+phcomplex.imDis

EndP
________________________________________________
Code: [Select]
Proc fft_calc:
    Arguments @nfactor, @InputX, @XFactor, @XtStruct, @step, @twids
    Local @kCounter, @SComplex, @MaxStep, @MaxNCount, @CurXtStruct, @CurXFactor, @CurSComplex, @CurXtStruct2
    Uses ecx, ebx, edi

    If D@nfactor = 1
        mov edi D@XFactor
        mov eax D@InputX
        fld R$eax | fstp R$edi+phcomplex.reDis
        fldz | fstp R$edi+phcomplex.imDis
        ExitP
    End_If

    mov ecx D@step | shl ecx 1 | mov D@MaxStep ecx
    mov eax D@nfactor | shr eax 1 | mov D@MaxNCount eax
    shl eax 4 | add eax D@XtStruct | mov D@SComplex eax

    call fft_calc D@MaxNCount, D@InputX, D@SComplex, D@XFactor, D@MaxStep, D@twids

    mov eax D@step | mov ecx D@InputX | lea ebx D$ecx+eax*8
    call fft_calc D@MaxNCount, ebx, D@XtStruct, D@XFactor, D@MaxStep, D@twids

    mov eax D@XtStruct | mov D@CurXtStruct eax | mov D@CurXtStruct2 eax
    mov eax D@XFactor | mov D@CurXFactor eax
    mov eax D@SComplex | mov D@CurSComplex eax

    mov D@kCounter 0
    mov ebx D@MaxNCount
    .While D@kCounter < ebx

        mov ecx D@kCounter | imul ecx D@Step | shl ecx 4 | add ecx D@twids
        call mult_complex D@CurXtStruct2, D@CurXtStruct, ecx
        call add_complex D@CurXFactor, D@CurSComplex, D@CurXtStruct2

        mov eax ebx | add eax D@kCounter | shl eax 4 | add eax D@XFactor
        call sub_complex eax, D@CurSComplex, D@CurXtStruct2

        add D@CurXtStruct Size_Of_phcomplex
        add D@CurXFactor Size_Of_phcomplex
        add D@CurSComplex Size_Of_phcomplex

        inc D@kCounter
    .End_While

EndP
_________________________________________________
Code: [Select]
Proc mult_complex:
    Arguments @Result, @AFactor, @BFActor
    Uses ebx, esi, edi

    mov ebx D@AFactor
    mov esi D@BFactor
    mov eax D@Result

    fld R$ebx+phcomplex.reDis | fmul R$esi+phcomplex.reDis
    fld R$ebx+phcomplex.imDis | fmul R$esi+phcomplex.imDis
    fsubp ST1 ST0
    fstp R$eax+phcomplex.reDis

    fld R$ebx+phcomplex.reDis | fmul R$esi+phcomplex.imDis
    fld R$ebx+phcomplex.imDis | fmul R$esi+phcomplex.reDis
    faddp ST1 ST0
    fstp R$eax+phcomplex.imDis

EndP
_________________________________________________

Code: [Select]
Proc add_complex:
    Arguments @Result, @AFactor, @BFActor
    Uses ebx, esi

    mov ebx D@AFactor
    mov esi D@BFactor
    mov eax D@Result

    fld R$ebx+phcomplex.reDis | fadd R$esi+phcomplex.reDis | fstp R$eax+phcomplex.reDis
    fld R$ebx+phcomplex.imDis | fadd R$esi+phcomplex.imDis | fstp R$eax+phcomplex.imDis

EndP
_________________________________________________

Code: [Select]
Proc sub_complex:
    Arguments @Result, @AFactor, @BFActor
    Uses ebx, esi

    mov ebx D@AFactor
    mov esi D@BFactor
    mov eax D@Result

    fld R$ebx+phcomplex.reDis | fsub R$esi+phcomplex.reDis | fstp R$eax+phcomplex.reDis
    fld R$ebx+phcomplex.imDis | fsub R$esi+phcomplex.imDis | fstp R$eax+phcomplex.imDis

EndP
_________________________________________________

Example of usage:

Code: [Select]
[Data1: R$ 1, 2, 3, 4, 5] ; Input

; Output
[Data2:
 Data2.Val1: R$ 0
 Data2.Val2: R$ 0
 Data2.Val3: R$ 0
 Data2.Val4: R$ 0
 Data2.Val5: R$ 0]

call fft Data1, 5, Data2




Coding in Assembly requires a mix of:
80% of brain, passion, intuition, creativity
10% of programming skills
10% of alcoholic levels in your blood.

My Code Sites:
http://rosasm.freeforums.org
http://winasm.tripod.com