Well, we can reconstruct A34 by doing the inverse transform on FFTA34.
Obviously this makes the test somewhat circular - though the fact that a
sensible-looking A34 (with integer valued elements) exists is encouraging...
OK, then, here goes for hopefully final versions of script (fftw.ijs), test
script (testfftw.ijs) and updated lab (fftw.ijt)
David
fftw.ijs:
NB. built from project: ~Addons/math/fftw/fftw
NB. z definitions:
script_z_ '~system\main\dll.ijs'
coclass 'jfftw'
fftw_z_=: (_1 & fftwnd_jfftw_) :. (1 & fftwnd_jfftw_)
ifftw_z_=: (1 & fftwnd_jfftw_) :. (_1 & fftwnd_jfftw_)
NB. fftw utils
NB.
NB. cd 15!:0
NB. clean clean numbers near 0
NB. info cover for wdinfo
NB. matchclean if clean x-y is all 0
zzero=: 1j1-1j1
3 : 0''
if. IFUNIX do.
DLL=: 'libfftw3.so.3'
else.
DLL=: '"',~'"',jpath '~addons\math\fftw\libfftw3-3.dll'
end.
)
FFTW_FORWARD=: _1
FFTW_BACKWARD=: 1
FFTW_ESTIMATE=: 6 (33 b.) 1
FFTW_MEASURE=: 0
FFTW_VERSION=: 3.2
cd=: 15!:0
info=: wdinfo @ ('FFTW'&;)
matchclean=: 0: *./ . = clean @ , @: -
NB. =========================================================
NB.*clean v clean numbers to tolerance (default 1e_10)
NB. sets values less than tolerance to 0
clean=: (1e_10&$:) : (4 : 0)
if. L. y do.
x clean each y
else.
if. 16 ~: 3!:0 y do.
y * x <: |y
else.
j./"1 y* x <: | y=. +.y
end.
end.
)
NB. fftw
NB. =========================================================
NB.*createplan v create a plan
NB. y = shape ; in ; out ; direction; flag
NB.
NB. direction = FFTW_FORWARD | FFTW_BACKWARD
NB. flag = FFTW_ESTIMATE | FFTW_MEASURE
createplan=: 3 : 0
'shape in out dir flag'=. y
assert dir e. FFTW_FORWARD,FFTW_BACKWARD
assert flag e. FFTW_ESTIMATE, FFTW_MEASURE
shape=. ,shape
cmd=. DLL,' fftw_plan_dft + x i *i *j *j i i'
0 pick cmd cd (#shape);shape;in;out;dir;flag
)
NB. =========================================================
NB.*destroyplan v destroy a plan
destroyplan=: 3 : 0
cmd=. DLL,' fftw_destroy_plan + n x'
1 [ cmd cd y
)
NB. =========================================================
NB.*fftwnd d n-dimensional FFT
NB. x = _1 forward
NB. 1 backward
NB. y = data
fftwnd=: 4 : 0
shp=. $y
if. 0 e. shp do. y return. end.
in=. zzero + , |: y
out=. in * 0
assert x e. _1 1
dir=. x
plan=. createplan shp;in;out;dir;FFTW_ESTIMATE
fftwexecute plan
destroyplan plan
res=. |: (|.shp) $ out
if. dir=1 do. res % */shp end.
)
NB. =========================================================
NB.*fftwexecute d one call to n-dimensional FFT
NB. y = plan
fftwexecute=: 3 : 0
cmd=. DLL,' fftw_execute + n x'
1 [ cmd cd y
)
testfftw.ijs
NB. test fftw system
NB.
NB. load this file with 0!:2
NB.
NB. 0!:2 <'\jx\addon\fftw\fftw\testfftw.ijs'
load '~addons\math\fftw\fftw.ijs'
clean=: 1e_10&$: : (4 : 0)
if. (3!:0 y) e. 16 16384 do.
j./"1 y* x <: | y=. +.y
else.
y * x <: |y
end.
)
dfft=: 3 : '+/ y * ^ (#y) %~ (- o. 0j2 ) * */~ i.#y'
matchclean=: 0: *./ . = clean @ , @: -
round=: [ * [: <. 0.5"_ + %~
A34=: ". ;._2 (0 : 0)
0 1 1j7 2j10
1 0j6 4j1 5j8
8j14 3j12 5j21 9j20
)
FFTA34=: ". ;._2 (0 : 0)
39j99 _10j_10 10.3661j18.8301 9.3205j_6.6795
_8j6 _33.0788j_17.2417 _13.9282j3.0526 8.634j10.1698
_1j13 _25.3205j_41.3205 24.0789j_39.7584 _0.0718j_35.0526
)
NB. =========================================================
NB. test dat = ifftw fftw dat
f=: matchclean if...@fftw
f 3
f i.8
f i.4096
f i.2 3 4
f ?.(6?.10)$100
NB. =========================================================
NB. test against direct calculation
f=: dfft matchclean fftw
f 3
f i.8
f i.132
NB. =========================================================
NB. known examples:
a=: 2 0 1 0 0 0 1 0
b=: 4 2 0 2 4 2 0 2
a matchclean ifftw b
b matchclean fftw a
FFTA34 matchclean 1e_4 round fftw A34
NB. =========================================================
NB. createplan, performance 'measure'
NB. In this case, must initialize input array _after_ creating plan.
NB. Note use of in-place modification to do this, and to re-use the plan
with
NB. new input.
in=: 12 $ 1j1 NB. Make sure elements are complex
values.
out=: in * 0
P=: createplan_jfftw_ 3 4;in;out;_1;0
in=: (,|:A34) a: } in
fftwexecute_jfftw_ P
FFTA34 matchclean 1e_4 round |: (4 3) $ out
in=: (2 * in) a: } in
fftwexecute_jfftw_ P
FFTA34 matchclean 1e_4 round -: |: (4 3) $ out
destroyplan_jfftw_ P
NB. createplan, performance 'estimate'
NB. Again note use of in-place modification when re-using plan with new
input.
in=: , |: A34
out=: in * 0
P=: createplan_jfftw_ 3 4;in;out;_1;64
fftwexecute_jfftw_ P
FFTA34 matchclean 1e_4 round |: (4 3) $ out
in=: (2 * in) a: } in
fftwexecute_jfftw_ P
FFTA34 matchclean 1e_4 round -: |: (4 3) $ out
destroyplan_jfftw_ P
fftw.ijt:
LABTITLE=: 'Fast Fourier Transform'
LABWINDOWED=: 0
LABFOCUS=: 0
NB. =========================================================
Lab Section FFTW
Fast Fourier Transform routines in J are based on the FFTW library, i.e. the
Fastest Fourier Transform in the West.
FFTW is a comprehensive collection of fast C routines for computing the
discrete Fourier transform (DFT) in one or more dimensions, of both real and
complex data, and of arbitrary input size.
FFTW was developed at MIT by Matteo Frigo and Steven G. Johnson, and a full
description is available at the FFTW home page:
http://theory.lcs.mit.edu/~fftw/
FFTW is licensed under the GNU General Public License, see
addons/math/fftw/gpl.htm.
)
NB. =========================================================
Lab Section Source Files
The FFTW routines are packaged in file libfftw3-3.dll, which includes the
complete set of complex number routines in FFTW.
These routines are accessible to J through the J DLL call mechanism as
documented in User Manual chapter "Dlls and Memory Management".
)
NB. =========================================================
Lab Section
The scripts and other files for FFTW are in subdirectory fftw/math
of the addons directory, i.e.
)
jpath '~addons/math/fftw'
NB. =========================================================
Lab Section Loading FFTW
Load the fftw.ijs script to access FFTW.
The main definitions for FFTW are read into locale jfftw.
Two verbs are defined in the z locale:
fftw forward FFT
ifftw backward FFT
)
load 'math/fftw'
NB. =========================================================
Lab Section Examples
Here is a simple example:
)
fftw i.8
NB. =========================================================
Lab Section
This example is easily verified by direct computation in J.
The discrete FFT for vector f is given by:
N-1
F(k) = sigma f(n) exp ((- j 2 pi n k)/N) k = 0 ... N-1
n=0
The verb dfft defined below implements the above expression, and matches the
result of fftw. However, it is very inefficient for large arguments.
Note that we remove small rounding errors by using the verb clean from the
numeric script.
)
dfft=: 3 : '+/ y * ^ (#y) %~ (- o. 0j2 ) * */~ i.#y'
require 'numeric'
clean (dfft - fftw) i.8
NB. =========================================================
Lab Section
The inverse FFT is given by:
ifftw
or equivalently:
fftw inverse
Note that the inverse is normalized, so that the inverse FFT returns the
original data.
)
ifftw fftw i.8
fftw inverse fftw i.8
NB. =========================================================
Lab Section
Doing a FFT followed by an inverse FFT returns the original data, so that
the following expression is everywhere zero:
dat - ifftw fftw dat
In practice, there may be small rounding errors. Again, we use clean to
remove these.
)
A=: ?.~ 8
A - ifftw fftw A
clean A - ifftw fftw A
NB. =========================================================
Lab Section
In general, the arguments to fftw and ifftw are complex, multidimensional
arrays.
The next section computes the FFT on a random complex, 6-dimensional array
of shape 3 4 5 6 7 8. This may take a few seconds to run.
)
NB. =========================================================
Lab Section
)
A=: j./ ?. 1 + i. 2 3 4 5 6 7 8
$A
$B=: fftw A
$C=: ifftw B
+/ , clean A - C NB. A matches C
NB. =========================================================
Lab Section Other Facilities
The verbs fftw and ifftw should cover practically all uses of FFTW.
However, there are two other facilities that are designed to optimize
repeated uses of the FFTW library on arguments of the same shape, which are
supported by the utilities in the jfftw locale:
saved plans
performance measures
)
NB. =========================================================
Lab Section
When using the underlying facilities directly, take care to ensure the
arguments given are correct. Calling DLL procedures incorrectly can crash
your system or corrupt memory.
)
NB. =========================================================
Lab Section Saved Plans
The verbs fftw and ifftw dynamically create and destroy plans. If you are
calling an FFT on same sized arrays repeatedly, you can create and reuse a
plan, destroying it when no longer required.
In the following example, a plan is created for a 2 3 4 array, and FFT is
called twice using it:
)
shape=: 2 3 4
data=: j./? 1+ i.2, shape NB. random data
NB. get data in right shape for library call and allocate
NB. space for output.
in=: (,|:data)
out=: in * 0
[P=: createplan_jfftw_ shape;in;out;_1;64 NB. create plan
NB. execute plan. Result goes to out, which we must reshape
fftwexecute_jfftw_ P
|: (|.shape) $ out
data=: j./? 1+ i.2, shape NB. new random data
in=: (,|:data) a: } in NB. in-place update of input
data
fftwexecute_jfftw_ P NB. reuse existing plan
|: (|.shape) $ out
destroyplan_jfftw_ P NB. destroy plan
NB. =========================================================
Lab Section Performance Measures
The final argument to the verb createplan is a flag, 0=MEASURE or
64=ESTIMATE, where ESTIMATE is the default.
Plans created with the MEASURE flag are optimized for a given array size. It
takes much longer to create such plans than with ESTIMATE, and so you would
do this only when a large number of FFTs are to be calculated for the same
array shape.
When creating plans with the ESTIMATE flag you may set up the input either
before or after creating the plan. But when creating plans with the MEASURE
flag, you must set up the input after creating the plan.
The next section creates a plan with a performance measure. It will take a
little longer to run, most of the time being spent creating the plan.
)
NB. =========================================================
Lab Section
)
shape=: 2 3 4
in=: (*/ shape) $ 1j1 NB. input must hold complex
values
out=: in * 0
[P=: createplan_jfftw_ shape;in;out;_1;0 NB. create plan
NB. populate input after creating plan.
data=: j./? 1+ i.2, shape
in=: (,|:data) a: } in
fftwexecute_jfftw_ P NB. execute plan
|: (|.shape) $ out NB. result goes to out
destroyplan_jfftw_ P NB. destroy plan
"bill lam" <[email protected]> wrote in message
news:[email protected]...
> lun, 22 Mar 2010, David Hotham skribis:
>> Hm. fftw A34 is still not equal to FFTA34.
>>
>> On win32 I see the same failure with the old script, though (and the new
>> script gives the same answer for fftw A34 as the old script did).
>>
>> Have you tried this yourself? Does the test pass for you (with either
>> old
>> or new script)?
>
> No, I didn't. I use 64-bit linux. It may employ some other methods to
> get the actual A34, such as run on some ancient J or ifftw2 or from
> matlab.
>
> --
> regards,
> ====================================================
> GPG key 1024D/4434BAB3 2008-08-24
> gpg --keyserver subkeys.pgp.net --recv-keys 4434BAB3
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
