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 :(
[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
_________________________________
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
________________________________________________
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
_________________________________________________
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
_________________________________________________
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
_________________________________________________
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:
[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