[futurebasic] Re: [FB] x-fb FFT and/or FIR

Message: < previous - next > : Reply : Subscribe : Cleanse
Home   : December 2001 : Group Archive : Group : All Groups

From: Ken Shmidheiser <k.shmidheiser@...>
Date: Wed, 5 Dec 2001 01:07:47 -0500
Concerning the ongoing thread about a port of the Fast Fourier 
Transform algorithm into FB^3, following its my translation of Murphy 
McCauley's Visual Basic FFT source code into FB^3. The original 
source code and discussion of FFT can be found here:

<http://www.intersrv.com/~dcross/fft.html#vb1>

Warning: This code is not tested and, in fact, I don't even know what 
it does. However, porting thousands of lines of VB code to FB^3 has 
given me a little understanding of the translation process, so this 
may be useful as springboard for further development.

Ken

Caution: Watch for e-mail line breaks and lost constant underscores 
on Associate server.

'-------- BEGIN FB^3 FAST FOURIER TRANSFORM CODE  ------

/*
    Based on:

    VB FFT Release 1.0
    by Murphy McCauley (MurphyMc@...)
    08/01/99

    This code is very, very heavily based on Don Cross's fourier.pas
    Turbo Pascal Unit for calculating the Fast Fourier Transform.
    I've not implemented all of his functions, though I may well do
    so in the future.

    FB^3 translation by Ken Shmidheiser, 12-4-01

*/

BEGIN GLOBALS
DIM DYNAMIC  RealIn( _maxInt ) AS DOUBLE
DIM DYNAMIC ImageIn( _maxInt ) AS DOUBLE
DIM DYNAMIC RealOut( _maxInt ) AS DOUBLE
DIM DYNAMIC ImagOut( _maxInt ) AS DOUBLE
END GLOBALS


CLEAR LOCAL FN NumberOfBitsNeeded( PowerOfTwo AS LONG )
DIM AS BYTE i

FOR i = 0 TO 16
LONG IF ( PowerOfTwo AND (2 ^ i) ) <> 0
FN NumberOfBitsNeeded = i
EXIT FN
END IF
NEXT i

END FN


CLEAR LOCAL FN IsPowerOfTwo( X AS LONG )
DIM AS BOOLEAN IsPowerOfTwo

IF (X < 2) THEN IsPowerOfTwo = _false: EXIT FN
IF (X AND (X - 1)) = _false THEN IsPowerOfTwo = _true

END FN = IsPowerOfTwo


CLEAR LOCAL FN ReverseBits( Index AS LONG, NumBits AS BYTE )
DIM AS LONG i, Rev

FOR i = 0 TO NumBits - 1
Rev = (Rev * 2) OR (Index AND 1)
Index = Index \ 2
NEXT i

END FN = Rev


CLEAR LOCAL FN Transform( NumSamples AS LONG )
DIM AS DOUBLE AngleNumerator, DeltaAngle, DeltaAr
DIM AS DOUBLE Alpha, Beta, TR, TI, AR, AI
DIM AS BYTE NumBits
DIM AS LONG I , j, k, n, BlockSize, BlockEnd

KILL DYNAMIC RealIn
KILL DYNAMIC ImageIn
KILL DYNAMIC RealOut
KILL DYNAMIC ImagOut

AngleNumerator = 2# * Pi

LONG IF ( FN IsPowerOfTwo( NumSamples ) = _false ) OR ( NumSamples < 2 )
PRINT "Error in procedure Fourier: NumSamples is "; NumSamples; ", 
which is not a positive integer power of two."
EXIT FN
END IF

NumBits = FN NumberOfBitsNeeded( NumSamples )
FOR I = 0 TO ( NumSamples - 1 )
j = FN ReverseBits( I, NumBits )
RealOut( j ) = RealIn( I )
ImagOut( j ) = ImageIn( I )
NEXT I

BlockEnd = 1
BlockSize = 2

WHILE BlockSize <= NumSamples
DeltaAngle = AngleNumerator / BlockSize
Alpha = SIN( 0.5 * DeltaAngle )
Alpha = 2# * Alpha * Alpha
Beta = SIN( DeltaAngle )

I = 0
WHILE I < NumSamples
AR = 1#
AI = 0#

j = I
FOR n = 0 TO BlockEnd - 1
k = j + BlockEnd
TR = AR * RealOut( k ) - AI * ImagOut( k )
TI = AI * RealOut( k ) + AR * ImagOut( k )
RealOut( k ) = RealOut( j ) - TR
ImagOut( k ) = ImagOut( j ) - TI
RealOut( j ) = RealOut( j ) + TR
ImagOut( j ) = ImagOut( j ) + TI
DeltaAr = Alpha * AR + Beta * AI
AI = AI - ( Alpha * AI - Beta * AR )
AR = AR - DeltaAr
j = j + 1
NEXT n

I = I + BlockSize
WEND

BlockEnd = BlockSize
BlockSize = BlockSize * 2
WEND

COMPRESS DYNAMIC RealIn
COMPRESS DYNAMIC ImageIn
COMPRESS DYNAMIC RealOut
COMPRESS DYNAMIC ImagOut

END FN

'----------- END FB^3 CODE -----------