- NICE.SOURCES (2:5030/1334.67) -------------------------------- NICE.SOURCES -
Msg : 22 из 707
From : Evgeniy Shvedun 2:463/283.211 Пон 18 Сен 00 15:55
To : All Срд 20 Сен 00 16:21
Subj : FFT
-------------------------------------------------------------------------------
Hi , All !
Мне нужно в потоке идущем на SB обpезать все выше 20KHz с очень хоpошим
качеством, для этого нужна пpавильная pеализация FFT
Как я понял мне надо:
1. Time domain to frequency domain
2. Обpезать
3. Frequency domain to time domain
Есть один pаботающий сыpец FFT на асме, но для 286:(
Кто-то с ним pаботал? Сомневаюсь в пpавильности pезультатов да и сpавнительно
тоpмозной он:(
Есть ли у All'а сыpец пpавильного и быстpого FFT на асме для 386? в идеале без
CoProcessor'а - очень скоpость кpитична.
Дайте плз очень нужно.
То-что есть:
=== Cut ===
; -------------------------------------
; Fast Fourier Transformation (FFT)
; -------------------------------------
;
; Filename : PCFFT.ASM
; Written by : J.G.G. Dobbe
; 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 (Turbo/Borland) C
; or (Turbo/Borlanc) Pascal.
; Comment : This module is written to link to:
;
; Turbo/Borland C
; Turbo/Borland Pascal
;
; 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!
;
; An assembler directive may be added to
; generate 80287 assembly code (instead of
; 8087 assembly code).
;
; TASM /dCOPROC287=1 PCFFT
;
; 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". The largest possible array is 16 K because an array
; variable may not axceed a segment boundary. The arrays may be
; placed in any arbitrary segment.
;
; 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)
;
; 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 segment
; named CODE. Local variables are also stored in segment CODE.
; ------------+- CoProcessor selection code -------------------------
IFDEF COPROC287
.286 ; For 80286 processor
.287 ; For 80287 coprocessor
ENDIF
; ------------+- 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
; Borland Pascal:
BdRe EQU DWORD PTR [BP+30] ; "Re"-array
BdIm EQU DWORD PTR [BP+26] ; "Im"-array
BdPower EQU [BP+24] ; "Pwr"
BdDir EQU [BP+22] ; "Dir"
ELSE
; Borland C:
BdRe EQU DWORD PTR [BP+22] ; "Re"-array
BdIm EQU DWORD PTR [BP+26] ; "Im"-array
BdPower EQU [BP+30] ; "Pwr"
BdDir EQU [BP+32] ; "Dir"
ENDIF
; ------------+- Start of code segment ------------------------------
CODE SEGMENT BYTE PUBLIC ; 'CODE' class name
ASSUME CS:CODE,DS:CODE ; ignored by TPascal
IFDEF PSCL ; Fft in asm
PUBLIC Fft
ELSE
PUBLIC _Fft
ENDIF
; ------------+- Variable definitions -------------------------------
AngleCounter DW (?) ; AngleCounter
FftLength DW (?) ; Length of FFT-arrs
FftPower DW (?) ; 2^FftPower = FftLen
FlyDistance DW (?) ; Fly-distance
Index1 DW (?) ; Index1
Index2 DW (?) ; Index2
Re1Pointer DW (?) ; Pnter to Re[Index1]
Re2Pointer DW (?) ; Pnter to Re[Index2]
Im1Pointer DW (?) ; Pnter to Im[Index1]
Im2Pointer DW (?) ; 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
;
Shuffle2Arr PROC FAR ;
PUSH AX ; CPU regs on stack
PUSH BX ;
PUSH CX ;
PUSH DX ;
PUSH SI ;
PUSH DI ;
PUSH ES ;
MOV DX,[FftPower] ; DL = Nr of bits shuffle-wrd
MOV DH,DL ; DH = Nr of bits shuffle-wrd
MOV CX,[FftLength] ; CX = counter
MOV SI,00H ; First index (IndexOld)
NextShuffleIndex: ;
MOV BX,SI ; Find shuffled word
MOV AX,00H ; BX = source index
CLC ; AX = destination index
NextB2: ; Clear Carry-flag
RCR BX,1 ; .-----------. .---.
RCL AX,1 ; BX=|..5 4 3 2 1|->-| C |--.
DEC DL ; `-----------' `---' |
JG NextB2 ; .-----------. |
MOV DL,DH ; AX=|1 2 3 4 5..|-<--------'
; `-----------'
CMP AX,SI ; IndexNew > IndexOld
JLE SkipIndex2 ; No? Skip index
PUSH DX ; Save DX
LES DI,BdRe ; Ptr to first element of a
MOV BX,SI ; Find IndexOld
SHL BX,1 ;
SHL BX,1 ;
PUSH BX ; Save for array b
ADD DI,BX ;
FLD DWORD PTR ES:[DI] ; 0 a[IndexOld]
LES BX,BdRe ; Ptr to first element of a
SHL AX,1 ; Find IndexNew
SHL AX,1 ;
PUSH AX ; Save for array b
ADD BX,AX ;
FLD DWORD PTR ES:[BX] ; 1 a[IndexNew]
FSTP DWORD PTR ES:[DI] ; a[IndexOld] = a[IndexNew]
FSTP DWORD PTR ES:[BX] ; a[IndexNew] = a[IndexOld]
POP AX ;
LES DI,BdIm ; Ptr to first element of b
POP BX ; Find IndexOld
ADD DI,BX ;
FLD DWORD PTR ES:[DI] ; 0 b[IndexOld]
LES BX,BdIm ; Ptr to first element of b
ADD BX,AX ;
FLD DWORD PTR ES:[BX] ; 1 b[IndexNew]
FSTP DWORD PTR ES:[DI] ; a[IndexOld] = a[IndexNew]
FSTP DWORD PTR ES:[BX] ; a[IndexNew] = a[IndexOld]
POP DX ; Restore DX
SkipIndex2:
INC SI ; Next index
LOOP NextShuffleIndex ;
POP ES ;
POP DI ;
POP SI ;
POP DX ;
POP CX ;
POP BX ;
POP AX ;
RET ; End of function Shuffle2Arr
Shuffle2Arr ENDP
; ------------+- NormalizeArr ---------------------------------------
;* * * * * * * * * * * * * * * * * * * * * * *
;* *
;* SUBROUTINE NormalizeArr(Var Arr:FftArr); *
;* *
;* * * * * * * * * * * * * * * * * * * * * * *
; Bitlength
; Normalizes array Arr by dividing each element by Sqrt(2 )
NormalizeArr PROC FAR
BdArr2 EQU DWORD PTR [BP+20] ; Borland C "Arr" pointer
PUSH AX ; Save CPU registers on stack
PUSH BX ;
PUSH CX ;
PUSH SI ;
PUSH DI ;
PUSH ES ;
PUSH DS ;
PUSH BP ;
MOV BP,SP ; BP = SP
MOV BX,[FftPower] ; Save BitLength for later
MOV CX,[FftLength] ;
LES SI,BdArr2 ; ES:SI points to Arr2
MOV AX,SEG ScaleFactor ; DS:DI points to scale-factor
MOV DS,AX ;
MOV DI,OFFSET ScaleFactor ;
DEC BX ;
SHL BX,1 ;
SHL BX,1 ;
ADD DI,BX ;
FLD DWORD PTR [DI] ; Scale-factor to coprocessor
ScaleElement:
FLD DWORD PTR ES:[SI] ; Array-element to coprocessor
FMUL ST,ST(1) ; ST(0) = ST(0) * ST(1)
FSTP DWORD PTR ES:[SI] ;
ADD SI,04H ; Next element
LOOP ScaleElement ; All elements had?
POP BP ; Restore CPU registers
POP DS ;
POP ES ;
POP DI ;
POP SI ;
POP CX ;
POP BX ;
POP AX ;
RET 04H ; End of function NormalizeArr
NormalizeArr ENDP
; ------------+- Fft ------------------------------------------------
;* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
;* *
;* void Fft(float *Re, float *Im, int Pwr, int Dir); *
;* *
;* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
;
; FlyCount CX
; section DX
; ------------+- Entry code -----------------------------------------
IFDEF PSCL
Fft PROC FAR ; Pascal entry code
ELSE
_Fft PROC FAR ; C entry code
ENDIF
; ------------+- Actual FFT algorithm -------------------------------
PUSH AX ; Save CPU registers on stack
PUSH BX ;
PUSH CX ;
PUSH DX ;
PUSH SI ;
PUSH DI ;
PUSH DS ;
PUSH ES ;
PUSH BP ;
MOV BP,SP ; BP = SP
MOV AX,CODE ; Load CODE segment
MOV DS,AX ;
FSAVE CoProcState ; Save CoProc status and init
MOV [CoCPUTemp],CoCPUCmd ; Initialize CoProcessor
FLDCW [CoCPUTemp] ;
MOV CX,BdPower ; Length array = 2
MOV [FftPower],CX ; Store power = WordLength
MOV AX,01H ;
power: ;
SHL AX,1 ;
LOOP power ;
MOV [FftLength],AX ;
CALL Shuffle2Arr ; Shuffle arrays Re and Im
MOV DX,01H ; Section = 1;
MOV AX,BdDir ; if(BdDir>=0)AngleCounter=1;
CMP AX,00H ; else AngleCounter=15;
JLE InverseFft ;
NormalFft:
MOV [AngleCounter],1 ; Normal FFT
JMP InitEnd ;
InverseFft:
MOV [AngleCounter],15 ; Inverse FFT
InitEnd:
MainLoop:
CMP DX,[FftLength] ; while (section < N)
JL SkipJump01 ;
JMP Normalize ;
SkipJump01:
MOV AX,DX ; FlyDistance=2 * section;
SHL AX,1 ;
MOV [FlyDistance],AX ;
MOV BX,OFFSET CosArray ; 0 c=CosArray[AngleCounter];
MOV AX,[AngleCounter] ;
DEC AX ;
SHL AX,1 ;
SHL AX,1 ;
ADD BX,AX ;
FLD DWORD PTR [BX] ;
MOV BX,OFFSET SinArray ; 1 s=SinArray[AngleCounter];
ADD BX,AX ;
FLD DWORD PTR [BX] ;
FLD1 ; 2 Qr=1.0;
FLDZ ; 3 Qi=0.0;
MOV CX,0FFFFH ; for (FlyCount=0;FlyCount <
; (section-1); FlyCount++)
ForLoop:
INC CX ;
MOV AX,DX ;
DEC AX ;
CMP CX,AX ;
JLE SkipJump03 ;
JMP ForLoopEnd ;
SkipJump03:
MOV [Index1],CX ; Index1=FlyCount;
Repeat1:
MOV AX,[Index1] ; Index2=Index1 + section;
MOV BX,AX ;
ADD AX,DX ;
MOV [Index2],AX ;
SHL BX,01H ; BX = offset to [Index1]
SHL BX,01H ;
SHL AX,01H ; AX = offset to [Index2]
SHL AX,01H ;
LES DI,BdIm ; Store pointer to Im[Index1]
MOV SI,DI ;
ADD DI,BX ;
MOV [Im1Pointer],DI ;
ADD SI,AX ; Store pointer to Im[Index2]
MOV [Im2Pointer],SI ;
FLD DWORD PTR ES:[SI] ; 4 Im[Index2] element->CoProc
FLD ST(0) ; 5 Im[index2]
FMUL ST(0),ST(2) ; 5 Im[Index2]*Qi
LES DI,BdRe ; Store pointer to Re[Index1]
MOV SI,DI ;
ADD DI,BX ;
MOV [Re1Pointer],DI ;
ADD SI,AX ; Store pointer to Re[Index2]
MOV [Re2Pointer],SI ;
FLD DWORD PTR ES:[SI] ; 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 *
;* * * * * * * * * * * * * *
LES SI,BdRe ; ES = segment of Re-array
MOV SI,[Re1Pointer] ;
FLD DWORD PTR [TempR] ; 0 TempR
FLD DWORD PTR ES:[SI] ; 1 Re[Index1]
FLD ST(0) ; 2 Re[Index1]
FSUB ST(0),ST(2) ; 2 Re[Index1] - TempR
MOV DI,[Re2Pointer] ;
FSTP DWORD PTR ES:[DI] ; Re[Index2]=Re[Index1]-TempR
;* * * * * * * * * * * * * *
;* Re1=Re[Index1] + TempR *
;* * * * * * * * * * * * * *
FADD ST(0),ST(1) ; 1 Re[Index1] + TempR
FSTP DWORD PTR ES:[SI] ; Re[Index1] += TempR
FSTP ST(0) ; Pop CoProcessor register
;* * * * * * * * * * * * * *
;* Im2=Im[Index1] - TempI *
;* * * * * * * * * * * * * *
LES SI,BdIm ; ES = segment of Im-array
MOV SI,[Im1Pointer] ;
FLD DWORD PTR [TempI] ; 0 TempI
FLD DWORD PTR ES:[SI] ; 1 Im[Index1]
FLD ST(0) ; 2 Im[Index1]
FSUB ST(0),ST(2) ; 2 Im[Index1] - TempI
MOV DI,[Im2Pointer] ;
FSTP DWORD PTR ES:[DI] ; Im[index2]=Im[Index1]-TempI
;* * * * * * * * * * * * * *
;* Im1=Im[Index1] + TempI *
;* * * * * * * * * * * * * *
FADD ST(0),ST(1) ; 1 Im[Index1] + TempI
FSTP DWORD PTR ES:[SI] ; Im[Index1] += TempI
FSTP ST(0) ; Pop CoProcessor register
MOV AX,[Index1] ; Index1 += FlyDistance;
MOV BX,[FlyDistance] ;
ADD AX,BX ;
MOV [Index1],AX ;
MOV BX,[FftLength] ; Index1 > (N-1)?
DEC BX ;
CMP AX,BX ;
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 DX,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 AX,BdDir ; if (BdDir == 1) NoNormalize;
CMP AX,01H ; else Normalize;
JE NoNormalize ;
LES SI,BdRe ; Normalize Re
PUSH ES ;
PUSH SI ;
CALL NormalizeArr ;
LES SI,BdIm ; Normalize Im
PUSH ES ;
PUSH SI ;
CALL NormalizeArr ;
NoNormalize:
FftEnd:
FRSTOR CoProcState ; Restore entire CoPU status
POP BP ; Restore CPU registers
POP ES ;
POP DS ;
POP DI ;
POP SI ;
POP DX ;
POP CX ;
POP BX ;
POP AX ;
; ------------+- 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 ------------------------------
CODE ENDS
END
=== Cut ===
With best wishes , Evgeniy.
... -= L.S.D. =- -= ZX SPECTRUM =- -= ASM =-
---
* Origin: Песок, опилки, бустилат - Ваш любимый Фрут энд Hат ! (2:463/283.211)
|
|
