Does results from fftw3 agree with that from fftw2? Can you also
verify fftw pass these tests?
NB. test fftw system
NB.
NB. load this file with 0!:2
NB.
NB. 0!:2 <'\jx\addon\fftw\fftw\testfftw.ijs'
load '~system\packages\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=: j./?.1+i.2 3 4
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 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 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
NB. existimate:
P=: createplan_jfftw_ 3 4;_1;0
FFTA34 matchclean 1e_4 round fftwndone_jfftw_ P;_1;A34
FFTA34 matchclean 1e_4 round -: fftwndone_jfftw_ P;_1;A34*2
destroyplan_jfftw_ P
NB. performance measure:
P=: createplan_jfftw_ 3 4;_1;1
FFTA34 matchclean 1e_4 round fftwndone_jfftw_ P;_1;A34
FFTA34 matchclean 1e_4 round -: fftwndone_jfftw_ P;_1;A34*2
destroyplan_jfftw_ P
dim, 21 Mar 2010, David Hotham skribis:
> "Therefore I suspect the problem originated from the window binary"
>
> It would be nice to put the blame elsewhere but, alas, this is not so.
>
> The issue seems to be that the definition of FFTW_ESTIMATE has changed
> between FFTW2 and FFTW3. What was zero is now (1U << 6).
>
> Here's a revised version of the script that:
> - incorporates your changes
> - removes some crud (or, at any rate, it seems to me to be crud; and
> removing it makes no difference on Windows)
> - fixes FFTW_ESTIMATE
>
> ... and (on w32 and w64) gives sensible results first time.
>
> David
>
>
> 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
> dir=. dir
> flag=. flag
> 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
> 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
> )
>
>
> "bill lam" <[email protected]> wrote in message
> news:[email protected]...
> [---=| TOFU protection by t-prot: 136 lines snipped |=---]
--
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
