[clean-list] 2 dimensional FFT code
Siegfried Gonzi
siegfried.gonzi@kfunigraz.ac.at
Thu, 08 Nov 2001 10:53:58 +0100
Enclosed a code for the two dimensional FFT (the one dimensional FFT is also posted again).
Timings:
This is strange. The following is valid for a 450MHz Pentium III machine with 128MB RAM.
A 1024x1240 2-dimensional FFT takes in the best cases 4sec. The problem is that the output results (the DOS window) are always slightly higher than the timing results (maybe the timing-profile overhead is taken into
account).
Yorick takes for a 1024x1204 complex transform 2sec. The Yorick times are *always* constant! With Clean it can happen that the program runs 10sec, or often 5sec; but 4sec are the best results.
The Yorick code itself is not very fast (compared to IDL or Matlab); okay IDL is very expensive I wouldn't be surprised that they have implemented a highly tuned one or even an assembler-version.
The Clean version could be 30% or even more faster, but due to the insane array types this will not happen. After the first transform I have to transpose the array (in order to access the columns effectively):
TransposeFFT:: !Int !*{#*{#Real}} -> !*{#*{#Real}}
For transposing I have to pass (or better: Clean expects) at least the innermost dimension as strict. BUT the first transform (before transposing in place):
fft_rows:: !Int !Int !Int !*{#{#Real}} !Real -> !*{#*{#Real}}
delivers the array ONLY in that format: !*{#*{#Real}}
when I use as exit strategy (in fft_rows):
| i == ( sizeA ) = { { el \\ el <-: array.[ij] } \\ ij <- [0..(dec sizeA )] }
BUT I WANT to use:
| i == (sizeA ) = array
NO, this is not allowed, because the compiler then puts back only an array of the following shape:
!*{#{#{Real}}
because this array format will be passsed to fft_rows and Clean is not willing to coarce this type to:
!*{#*{#Real}}
but this type is used in the transpose function.
I tried to work-around the problem, but IT IS IMPOSSIBLE.
And this costs Clean 30% and even more.
S. Gonzi
==
module fft_save
import StdEnv
/////////////////////////////////////////////////////////////////////////////////
//The program comes with no warranty!
//
//The one dimensional FFT is based on the Numerical Recipes
// one. To call a one dimensional FFT, you have to pass
// a complex array, where the first value
//holds the real part, the second value the imaginary part
//the third the real, and so forth.
//
//e.g. Start= FFT1d 8 (BitRev 4 {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0}) 1.0
//
//will transform a complex array of the the following shape:
//
// {0.0+1.0i,2.0+3.0i,4.0+5.0i,6.0+7.0i}
//
//the dimension here is 8 (and not 4)
//
//1.0 is the sign for the forward 1.0 or backward (-1.0) transform
//
// Before you pass the array to FFT1d you have to make the bit-reversing
//part: BitRev 8 array
//
//where 8 denotes still the dimension
//
// THE results are not scaled for N; if you nedd absolute values you may
//have the pleasure to do it yourself and divide the output by the dimension
//
//(C) Siegfied Gonzi 2001
/////////////////////////////////////////////////////////////////////////////////////////////////////
FFT1d:: !Int !*{#Real} !Real -> !*{#Real}
FFT1d n marray isign = LOOP1 2 marray
where
LOOP1:: !Int !*{#Real} -> !*{#Real}
LOOP1 mmax marray
| mmax>=n = marray
= LOOP1 ( mmax*2 ) ( LOOP2 1 marray 1.0 0.0 )
where
istep = mmax*2
theta = isign*( 6.2831853071959 / toReal( mmax ) )
wtemp = sin(0.5*theta)
wpr = ( -2.0 )*wtemp*wtemp
wpi = sin( theta )
LOOP2:: !Int !*{#Real} !Real !Real -> !*{#Real}
LOOP2 m marray wr wi
| m>=( mmax ) = marray
#! wtemp = wr
= LOOP2 ( m+2 ) ( LOOP3 m marray ) ( wr+wtemp*wpr-wi*wpi ) ( wi+wi*wpr+wtemp*wpi )
where
LOOP3:: !Int !*{#Real} -> !*{#Real}
LOOP3 i marray
| i>=n = marray
#! xi = marray.[i-1]
#! yi = marray.[i]
#! xj = marray.[i-1+mmax]
#! yj = marray.[i+mmax]
#! tempr = wr*xj - wi*yj
#! tempi = wr*yj + wi*xj
= LOOP3 ( i+istep ) { marray &
[i-1+mmax] = ( xi - tempr ),
[i+mmax] = ( yi - tempi ),
[i-1] = ( xi + tempr ),
[i] = ( yi + tempi )
}
/////////////////////////////////////////////////////////////////////////////
//The bit reversing part
//The first number denotes the dimension (imaginary
// and real)!
//This program has to be used before the FFT1d
///////////////////////////////////////////////////////////////////////////
BitRev:: !Int !*{#Real} -> !*{#Real}
BitRev n marray = BitShuffling 0 0 marray
where
BitShuffling:: !Int !Int !*{#Real} -> !*{#Real}
BitShuffling i j marray
| i>=( n-1 ) = marray
| j>i
#! xj = marray.[j]
#! xi =marray.[i]
#! yi = marray.[i+1]
#! yj = marray.[j+1]
= BitShuffling ( i+2 ) ( m1 ) { marray & [i] = xj, [j] = xi, [j+1] = yi, [i+1] = yj }
| otherwise = BitShuffling ( i+2 ) ( m1 ) marray
where
m1 = ( LOOP2 ( n/2 ) j )
LOOP2:: !Int !Int -> !Int
LOOP2 m j
| m <2 || j<m = ( j+m )
= LOOP2 ( m/2 ) ( j-m )
//////////////////////////////////////////////////////////////////////////////////////////////
//The 2-dimensional FFT is not based on any code (it is
//actually a self-creation by Gonzi); but it relies on the FFT1d
//
//The array to pass has the following shape:
//
// {{1.0,2.0i,3.0,4.0i...first row with real and complex numbers},
// {1.0,2.0i,3.0,4.0i,...second row with real and complex numbers},
// {last row with real and complex numbers}}
//
//The output orientatiom is as the input one. In all cases the number of
// columns is *two times* that of rows; because every columns
//holds real *and* imaginary parts. THIS is not checked for! You
//have to take care.
//
//THE Results are not scaled for the dimension (like Yorick) ;
//if you need it, you can divide the output by the dimension
//
//(C) Siegfried Gonzi 2001
/////////////////////////////////////////////////////////////////////////////////////////////////////
FFT2d:: !*{#*{#Real}} !Real -> !*{#{#Real}}
FFT2d marray sign
#! sizeA = size marray
#! sizeR = size marray.[0]
#! marray_row = ( fft_rows 0 sizeA sizeR marray sign )
#! trans_array = ( TransposeFFT sizeA marray_row )
#! marray_col = ( fft_rows 0 sizeA sizeR trans_array sign )
= ( TransposeFFT sizeA marray_col )
where
fft_rows:: !Int !Int !Int !*{#{#Real}} !Real -> !*{#*{#Real}}
fft_rows i sizeA sizeR array sign
| i == ( sizeA ) = { { el \\ el <-: array.[ij] } \\ ij <- [0..(dec sizeA )] }
//| i == (sizeA ) = array
#! row = array.[i]
#! bitRev = ( BitRev sizeR { el \\ el <-: row } )
#! fft1d = ( FFT1d sizeR bitRev 1.0 )
= fft_rows ( inc i ) sizeA sizeR { array & [i] = fft1d } sign
//
TransposeFFT:: !Int !*{#*{#Real}} -> !*{#*{#Real}}
TransposeFFT sizeA array = Transpose_j 0 array
where
Transpose_j:: !Int !*{#*{#Real}} -> !*{#*{#Real}}
Transpose_j i array
| i == ( sizeA -1 ) = array
#! trans_col = ( Transpose_i (inc i) (2*i) i array )
= Transpose_j ( inc i ) trans_col
//
Transpose_i:: !Int !Int !Int !*{#*{#Real}} -> !*{#*{#Real}}
Transpose_i i j ij array
| i == ( sizeA ) = array
#! (a1,array) = array![i,j]
#! (a2,array) = array![i,(inc j)]
#! (a3,array) = array![ij,(2*i)]
#! (a4,array) = array![ij,(2*i+1)]
#! ar = swapBug a1 a2 ij i array
= Transpose_i ( inc i ) j ij { ar & [i,j] = a3, [i,(inc j)] = a4 }
where
swapBug:: !Real !Real !Int !Int !*{#*{#Real}} -> !*{#*{#Real}}
swapBug val1 val2 jc i array = { array & [jc,(2*i)] = val1, [jc,(2*i+1)] = val2 }
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//TEST
//
//A complex array of the shape:
//
// 0.0+1.0i,2.0+3.0i,4.0+5.i,6.0+7.0i
// 0.0+1.0i,2.0+3.0i,4.0+5.i,6.0+7.0
// 0.0+1.0i,2.0+3.0i,4.0+5.i,6.0+7.0
// 0.0+1.0i,2.0+3.0i,4.0+5.i,6.0+7.0
//
Start = FFT2d ( MyArray n1 n2 ) ( 1.0 )
where
n1 = 8
n2 = 4
//
MyArray:: !Int !Int -> !*{#*{#Real}}
MyArray n1 n2 = { { toReal(x) \\ x <- [0..(n1-1)] } \\ y <- [0..(n2-1)] }
//Yorick delivers here:
//
// [[48+64i,-32+0i,-16-16i,0-32i],[0+0i,0+0i,0+0i,0+0i],[0+0i,0+0i,0+0i,0+0i],[0+0i,0+0i,0+0i,0+0i]]
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// A complex 1024x1024 array
Start:: Real
Start = TakeOne( FFT2d ( MyArray n1 n2 ) ( -1.0 ) )
where
n1 = 2048
n2 = 1024
//
TakeOne:: !{#{#Real}} -> !Real
TakeOne array = array.[0,0]
//
MyArray:: !Int !Int -> !*{#*{#Real}}
MyArray n1 n2 = { { toReal(x) \\ x <- [0..(n1-1)] } \\ y <- [0..(n2-1)] }
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
==