home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
OS/2 Shareware BBS: 10 Tools
/
10-Tools.zip
/
fft32_a.zip
/
PCFFT32.ASM
< prev
next >
Wrap
Assembly Source File
|
1995-01-20
|
25KB
|
590 lines
.386
.387
.model flat
; -------------------------------------
; Fast Fourier Transformation (FFT)
; -------------------------------------
;
; Filename : PCFFT.ASM
; Written by : J.G.G. Dobbe
; : 32 bit adaptation by R. Canup
; Function : Fast Fourier Transform S/W. Two arrays hol-
; ding short-real floating point numbers (4
; bytes per real value) are used to perform
; the FFT after a call from C or Pascal compilers.
; Comment : This module is written to link to:
;
; Gnu C 2.5.4. 32 bit OS/2 compiler 1/19/95
;
; An assembler switch (/dPSCL) must be entered
; if the object must be used with Pascal. It
; generates the required Pacal calling con-
; vention.
;
; Note: A coprocessor must be present when
; using this code!
;
; The C prototype for this routine is:
;
; void Fft(float *Re, float *Im, int Pwr, int Dir);
;
; Performs Fast Fourier Transform on the arrays Re and Im. The
; length of the arrays is given as a power of two in variable
; "Pwr". Example: if array has 1024 values "Pwr" is 10 - since two to the
; 10th power is 1024.
;
; Variable "Dir" determines whether a normal or an inversed FFT
; will be performed:
;
; Dir >= 1 => FFT (Time domain to frequency domain)
; Dir <= 0 => Inversed FFT (Frequency domain to time domain)
;
; Example:
;
; Fft( &Re, &Im, 12, -1);
;
; This would do an inverse fft on a pair of arrays with 4096 floating point
; values.
;
; No error-checking is performed on overflow. You can be sure no
; overflow occurs when each element of Re and Im, is smaller than
; the maximum value of:
;
; 3.4 E 38 / N (N = Length of arrays Re and Im)
;
; The executable code for function Fft is placed in a common code segment.
; Local variables are stored in a data segment. Despite the
; propaganda that programs in OS/2 are all in one segment - the truth is that
; the Intel architecture won't allow writing to a code segment. In actuality
; the OS/2 model is more like a giant "small" model in a dos compiler: one
; segment for code, and one for data and stack. Since each of these segments
; has a limit of 1bffffffh (When a program bombs examine the data registers and
; segment limits) which is 448 megabytes - this is not much of a constraint. I
; would be very interested in seeing a program which coudn't fit in 448 mega-
; bytes.
; OS/2 traps out to an error if the calculations produce a denormalized result
; with the following initialization values.
; -------------- Constants used by FFT-main routine -----------------
; CoProcessor initialization:
;
; Bit Code Use
; 0 IM Invalid operation exception mask 0 Enabled
; 1 DM Denormalized operand exception mask 0 Enabled
; 2 ZM Zerodivide exception mask 0 Enabled
; 3 OM Overflow exception mask 0 Enabled
; 4 UM Underflow exception mask 1 Disabled
; 5 PM Precision exception mask 1 Disabled
; 6 - - - -
; 7 IEM Interrupt enable mask 0 Enabled
; 8 PC Precision control\ 1 - 64 bits
; 9 PC Precision control/ 1 /
; 10 RC Rounding control\ 0 - Nearest even
; 11 RC Rounding control/ 0 /
; 12 IC Infinity control 1 Affine
; 13 - - - -
; 14 - - - -
; 15 - - - -
CoCPUCmd EQU 1330H ; CoCPU init code
IFDEF PSCL
; These values may need changed if they are used with a 32 bit OS/2 Pascal
; Borland Pascal for dos values:
BdRe EQU DWORD PTR [EBP+30] ; "Re"-array
BdIm EQU DWORD PTR [EBP+26] ; "Im"-array
BdPower EQU [EBP+24] ; "Pwr"
BdDir EQU [EBP+22] ; "Dir"
ELSE
; The stack frame is as follows for Gnu C 2.5.4. :
;stack_frame STRUC
; redi dd ? ; register set caused by PUSHAD instruction
; resi dd ? ; .
; rebp dd ? ; .
; resp dd ? ; .
; rebx dd ? ; .
; redx dd ? ; .
; recx dd ? ; .
; reax dd ? ; .
; retaddr dd ? ; return address
; BdRe dd ? ; "Re" - array address
; BdIm dd ? ; "Im" - array address
; BdPower dd ? ; "Pwr"- 32 bit integer value
; BdDir dd ? ; "Dir"- 32 bit integer value
;stack_frame ENDS
;
; Thus the offsets for the values are:
;
BdRe EQU DWORD PTR [EBP+36] ; "Re"-array
BdIm EQU DWORD PTR [EBP+40] ; "Im"-array
BdPower EQU DWORD PTR [EBP+44] ; "Pwr"
BdDir EQU DWORD PTR [EBP+48] ; "Dir"
;
; These defines allow constructs such as:
; MOV EAX,BdRe
; instead of
; MOV EAX,DWORD PTR [EBP+36]
ENDIF
.data
AngleCounter DD (?) ; AngleCounter
FftLength DD (?) ; Length of FFT-arrs
FftPower DD (?) ; 2^FftPower = FftLen
FlyDistance DD (?) ; Fly-distance
Index1 DD (?) ; Index1
Index2 DD (?) ; Index2
Re1Pointer DD (?) ; Pnter to Re[Index1]
Re2Pointer DD (?) ; Pnter to Re[Index2]
Im1Pointer DD (?) ; Pnter to Im[Index1]
Im2Pointer DD (?) ; Pnter to Im[Index2]
TempR DD (?) ; TempR
TempI DD (?) ; TempI
CoCPUTemp DW (?) ; Temp data for CoCPU
CoProcState DB 94 DUP (?) ; CoProc status area
CosArray DD -1.00000000000000 ; TurboPascal doesn't
DD 0.00000000000000 ; accept this initia-
DD 0.70710678118655 ; lized data to be in
DD 0.92387953251129 ; the DATA segment
DD 0.98078528040323
DD 0.99518472667220
DD 0.99879545620517
DD 0.99969881869620
DD 0.99992470183914
DD 0.99998117528260
DD 0.99999529380958
DD 0.99999882345170
DD 0.99999970586288
DD 0.99999992646572
DD -1.00000000000000
DD 0.00000000000000
DD 0.70710678118655
DD 0.92387953251129
DD 0.98078528040323
DD 0.99518472667220
DD 0.99879545620517
DD 0.99969881869620
DD 0.99992470183914
DD 0.99998117528260
DD 0.99999529380958
DD 0.99999882345170
DD 0.99999970586288
DD 0.99999992646572
SinArray DD 0.00000000000000
DD -1.00000000000000
DD -0.70710678118655
DD -0.38268343236509
DD -0.19509032201613
DD -0.09801714032956
DD -0.04906767432742
DD -0.02454122852291
DD -0.01227153828572
DD -0.00613588464915
DD -0.00306795676297
DD -0.00153398018628
DD -0.00076699031874
DD -0.00038349518757
DD 0.00000000000000
DD 1.00000000000000
DD 0.70710678118655
DD 0.38268343236509
DD 0.19509032201613
DD 0.09801714032956
DD 0.04906767432742
DD 0.02454122852291
DD 0.01227153828572
DD 0.00613588464915
DD 0.00306795676297
DD 0.00153398018628
DD 0.00076699031874
DD 0.00038349518757
ScaleFactor DD 0.50000000000000
DD 0.25000000000000
DD 0.12500000000000
DD 0.06250000000000
DD 0.03125000000000
DD 0.01562500000000
DD 0.00781250000000
DD 0.00390625000000
DD 0.00195312500000
DD 0.00097656250000
DD 0.00048828125000
DD 0.00024414062500
DD 0.00012207031250
DD 0.00006103515625
; -------------- Shuffle2Arr ----------------------------------------
;* * * * * * * * * * * * * *
;* *
;* SUBROUTINE Shuffle2Arr; *
;* *
;* * * * * * * * * * * * * *
;
; Shuffles array a and b:
;
; index = binair --> shuffled index binair index
;
; 0 000 000 0
; 1 001 100 4
; 2 010 010 2
; 3 011 110 6
; 4 100 001 1
; 5 101 101 5
; 6 110 011 3
; 7 111 111 7
;
; 2
; n = word-length (3 in this example) => n = log (Length Array a)
; n must be > 0
;
; -------------- Start of code segment ------------------------------
.code
IFDEF PSCL ; Fft in asm
PUBLIC Fft
ELSE
PUBLIC _Fft
ENDIF
; -------------- Variable definitions -------------------------------
Shuffle2Arr PROC NEAR
PUSH EAX ; CPU regs on stack
PUSH EBX ;
PUSH ECX ;
PUSH EDX ;
PUSH ESI ;
PUSH EDI ;
; IN 32 bits we don't touch segments
MOV EDX,[FftPower] ; DL = Nr of bits shuffle-wrd
MOV DH,DL ; DH = Nr of bits shuffle-wrd
MOV ECX,[FftLength] ; ECX = counter
MOV ESI,00H ; First index (IndexOld)
NextShuffleIndex:
MOV EBX,ESI ; Find shuffled word
MOV EAX,00H ; EBX = source index
CLC ; EAX = destination index
NextB2:
RCR EBX,1 ; ,-----------, ,---,
RCL EAX,1 ;EBX=|,,5 4 3 2 1|->-| C |--,
DEC DL ; `-----------' `---' |
JG NextB2 ; ,-----------, |
MOV DL,DH ;EAX=|1 2 3 4 5,,|-<--------'
; `-----------'
CMP EAX,ESI ; IndexNew > IndexOld
JLE SkipIndex2 ; No? Skip index
PUSH EDX ; SAVE EDX
MOV EDI,BdRe ; Ptr to first element of a
MOV EBX,ESI ; Find IndexOld
SHL EBX,2 ;
PUSH EBX ; Save for array b
ADD EDI,EBX ;
FLD DWORD PTR [EDI] ; 0 a[IndexOld]
MOV EBX,BdRe ; Ptr to first element of a
SHL EAX,2 ; Find IndexNew
PUSH EAX ; Save for array b
ADD EBX,EAX ;
FLD DWORD PTR [EBX] ; 1 a[IndexNew]
FSTP DWORD PTR [EDI] ; a[IndexOld] = a[IndexNew]
FSTP DWORD PTR [EBX] ; a[Indexnew] = a[IndexOld]
POP EAX ;
MOV EDI,BdIm ; Ptr to first element of b
POP EBX ; Find IndexOld
ADD EDI,EBX ;
FLD DWORD PTR [EDI] ; 0 b[IndexOld]
MOV EBX,BdIm ; Prt to first element of b
ADD EBX,EAX ;
FLD DWORD PTR [EBX] ; 1 b[IndexNew]
FSTP DWORD PTR [EDI] ; a[IndexOld] = a[IndexNew]
FSTP DWORD PTR [EBX] ; a[IndexNew] = a[IndexOld]
POP EDX ; Restore EDX
SkipIndex2:
INC ESI ; Next index
LOOP NextShuffleIndex ;
POP EDI ;
POP ESI ;
POP EDX ;
POP ECX ;
POP EBX ;
POP EAX ;
RET ; End of function Shuffle2Arr
Shuffle2Arr ENDP
; --------------- NormalizeArr --------------------------------------
;* * * * * * * * * * * * * * * * * * * * * * *
;* *
;* SUBROUTINE NomalizeArr(Var Arr:FftArr); *
;* *
;* * * * * * * * * * * * * * * * * * * * * * *
; Bitlength
; Normalizes array Arr by dividing each element by Sqrt(2 )
NormalizeArr PROC NEAR
BdArr2 EQU DWORD PTR [EBP+36] ; GNU C 2.5.4. "Arr" pointer
PUSHAD ; SAVE ALL THE REGISTERS
MOV EBP,ESP ; SET UP BASE POINTER
MOV EBX,[FftPower] ; Save BitLength for later
MOV ECX,[FftLength] ;
MOV ESI,BdArr2 ; ESI points to Arr2
MOV EDI,OFFSET ScaleFactor ;
DEC EBX ;
SHL EBX,2 ; Take advantage of barrel shifter
ADD EDI,EBX ;
FLD DWORD PTR [EDI] ; Scale-factor to coprocessor
ScaleElement:
FLD DWORD PTR [ESI] ; Array-element to coprocessor
FMUL ST,ST(1) ; ST(0) = ST(0) * ST(1)
FSTP DWORD PTR [ESI] ;
ADD ESI,04H ; Next element
LOOP ScaleElement ; All elemnts had?
POPAD ; Restore the registers
RET 04H ; End of function NormalizeArr
NormalizeArr ENDP
; -------------- Fft ------------------------------------------------
;* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
;* *
;* void Fft(float *Re, float *Im, int Pwr, int Dir); *
;* *
;* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
;
; FlyCount ECX
; section EDX
; -------------- Entry code -----------------------------------------
IFDEF PSCL
Fft PROC ; Pascal entry code
ELSE
_Fft PROC ; C entry code
ENDIF
; -------------- Actual FFT algorithm -------------------------------
PUSHAD ; Save all the registers on the stack
MOV EBP,ESP ; EBP = ESP
FSAVE CoProcState ; Save CoProc status and init
MOV [CoCPUTemp],CoCPUCmd ; Initialize CoProcessor
FLDCW [CoCPUTemp] ;
MOV ECX,BdPower ; Length array = 2
MOV [FftPower],ECX ; Store power = WordLength
MOV EAX,01H ;
power: ; ecx
SHL EAX,CL ; Create actual size of array 2
MOV [FftLength],EAX ;
CALL Shuffle2Arr ; Shuffle arrays Re and Im
MOV EDX,01H ; Section = 1;
MOV EAX,BdDir ; if(BdDir>=0)AngleCounter=1;
CMP EAX,00H ; else AngleCounter=15;
JLE InverseFft ;
NormalFft:
MOV [AngleCounter],1 ; Normal FFT
JMP InitEnd ;
InverseFft:
MOV [AngleCounter],15 ; Inverse FFT
InitEnd:
MainLoop:
CMP EDX,[FftLength] ; while (section < N)
JL SkipJump01 ;
JMP Normalize ;
SkipJump01:
MOV EAX,EDX ; FlyDistance=2 * section;
SHL EAX,1 ;
MOV [FlyDistance],EAX ;
MOV EBX,OFFSET CosArray ; 0 c=CosArray[AngleCounter];
MOV EAX,[AngleCounter] ;
DEC EAX ;
SHL EAX,2 ;
ADD EBX,EAX ;
FLD DWORD PTR [EBX] ;
MOV EBX,OFFSET SinArray ; 1 s=SinArray[AngleCounter];
ADD EBX,EAX ;
FLD DWORD PTR [EBX] ;
FLD1 ; 2 Qr=1.0;
FLDZ ; 3 Qi=0.0;
MOV ECX,0FFFFFFFFH ; for (FlyCount=0;FlyCount <
; (section-1); FlyCount++)
ForLoop:
INC ECX ;
MOV EAX,EDX ;
DEC EAX ;
CMP ECX,EAX ;
JLE SkipJump03 ;
JMP ForLoopEnd ;
SkipJump03:
MOV [Index1],ECX ; Index1=FlyCount;
Repeat1:
MOV EAX,[Index1] ; Index2=Index1 + section;
MOV EBX,EAX ;
ADD EAX,EDX ;
MOV [Index2],EAX ;
SHL EBX,02H ; EBX = offset to [Index1]
SHL EAX,02H ;
MOV EDI,BdIm ; Store pointer to Im[Index1]
MOV ESI,EDI ;
ADD EDI,EBX ;
MOV [Im1Pointer],EDI ;
ADD ESI,EAX ; Store pointer to Im[Index2]
MOV [Im2Pointer],ESI ;
FLD DWORD PTR [ESI] ; 4 Im[Index2] element->CoProc
FLD ST(0) ; 5 Im[index2]
FMUL ST(0),ST(2) ; 5 Im[Index2]*Qi
MOV EDI,BdRe ; Store pointer to Re[Index1]
MOV ESI,EDI ;
ADD EDI,EBX ;
MOV [Re1Pointer],EDI ;
ADD ESI,EAX ; Store pointer to Re[Index2]
MOV [Re2Pointer],ESI ;
FLD DWORD PTR [ESI] ; 6 Re[Index2]
FLD ST(0) ; 7 Re[Index2]
FMUL ST(0),ST(5) ; 7 Re[Index2]*Qr
FSUB ST(0),ST(2) ; 7 Re[Index2]*Qr-Im[Index2]*Qi
FSTP DWORD PTR TempR ; Store at TempR and pop
FMUL ST(0),ST(3) ; 6 Re[Index2]*Qi
FLD ST(2) ; 7 Im[Index2]
FMUL ST(0),ST(5) ; 7 Im[Index2]*Qr
FADD ST(0),ST(1) ; 7 Im[Index2]*Qr+Re[Index2]*Qi
FSTP DWORD PTR TempI ; Store at TempI and pop
FSTP ST(0) ; Pop CoProcessor registers
FSTP ST(0) ;
FSTP ST(0) ;
;* * * * * * * * * * * * * *
;* Re2=Re[Index1] - TempR *
;* * * * * * * * * * * * * *
MOV ESI,BdRe ; Re-array
MOV ESI,[Re1Pointer] ;
FLD DWORD PTR [TempR] ; 0 TempR
FLD DWORD PTR [ESI] ; 1 Re[Index1]
FLD ST(0) ; 2 Re[Index1]
FSUB ST(0),ST(2) ; 2 Re[Index1] - TempR
MOV EDI,[Re2Pointer] ;
FSTP DWORD PTR [EDI] ; Re[Index2]=Re[Index1]-TempR
;* * * * * * * * * * * * * *
;* Re1=Re[Index1] + TempR *
;* * * * * * * * * * * * * *
FADD ST(0),ST(1) ; 1 Re[Index1] + TempR
FSTP DWORD PTR [ESI] ; Re[Index1] += TempR
FSTP ST(0) ; Pop CoProcessor register
;* * * * * * * * * * * * * *
;* Im2=Im[Index1] - TempI *
;* * * * * * * * * * * * * *
MOV ESI,BdIm ; ES = segment of Im-array
MOV ESI,[Im1Pointer] ;
FLD DWORD PTR [TempI] ; 0 TempI
FLD DWORD PTR [ESI] ; 1 Im[Index1]
FLD ST(0) ; 2 Im[Index1]
FSUB ST(0),ST(2) ; 2 Im[Index1] - TempI
MOV EDI,[Im2Pointer] ;
FSTP DWORD PTR [EDI] ; Im[index2]=Im[Index1]-TempI
;* * * * * * * * * * * * * *
;* Im1=Im[Index1] + TempI *
;* * * * * * * * * * * * * *
FADD ST(0),ST(1) ; 1 Im[Index1] + TempI
FSTP DWORD PTR [ESI] ; Im[Index1] += TempI
FSTP ST(0) ; Pop CoProcessor register
MOV EAX,[Index1] ; Index1 += FlyDistance;
MOV EBX,[FlyDistance] ;
ADD EAX,EBX ;
MOV [Index1],EAX ;
MOV EBX,[FftLength] ; Index1 > (N-1)?
DEC EBX ;
CMP EAX,EBX ;
JG SkipJump02 ; No? To Repeat1
JMP Repeat1 ;
SkipJump02:
FLD ST(0) ; 4 Qi
FMUL ST(0),ST(3) ; 4 Qi*s
FLD ST(2) ; 5 Qr
FMUL ST(0),ST(4) ; 5 Qr*s
FLD ST(3) ; 6 Qr
FMUL ST(0),ST(6) ; 6 Qr*c
FSUB ST(0),ST(2) ; 6 Qr*c - Qi*s
FSTP ST(4) ; Store and pop
FLD ST(2) ; 6 Qi
FMUL ST(0),ST(6) ; 6 Qi*c
FADD ST(0),ST(1) ; 6 Qi*c + Qr*s
FSTP ST(3) ; store and pop
FSTP ST(0) ; Pop CoProcessor registers
FSTP ST(0) ;
JMP ForLoop ;
ForLoopEnd:
SHL EDX,01H ; section = section * 2;
INC [AngleCounter] ; AngleCounter=AngleCounter+1;
FSTP ST(0) ; Pop CoProcessor registers
FSTP ST(0) ;
FSTP ST(0) ;
FSTP ST(0) ;
JMP MainLoop ;
Normalize:
MOV EAX,BdDir ; if (BdDir == 1) NoNormalize;
CMP EAX,01H ; else Normalize;
JE NoNormalize ;
MOV ESI,BdRe ; Normalize Re
PUSH ESI ;
CALL NormalizeArr ;
MOV ESI,BdIm ; Normalize Im
PUSH ESI ;
CALL NormalizeArr ;
NoNormalize:
FftEnd:
FRSTOR CoProcState ; Restore entire CoPU status
POPAD ; Restore CPU registers
; -------------- Termination code -----------------------------------
IFDEF PSCL ; Termination code for Pascal
RET 0CH ; End of function Fft
Fft ENDP
ELSE ; Termination code for C
RET ; End of function Fft
_Fft ENDP
ENDIF
; -------------- Assembler termination ------------------------------
END
