Running on an ancient version.
9!:14''
4.06/2001-05-22/23:00
]A34=: j./?.1+i.2 3 4
0 1 1j7 2j10
1 0j6 4j1 5j8
8j14 3j12 5j21 9j20
The assumption on A34 should be correct. Thank you for your work!
mar, 23 Mar 2010, David Hotham skribis:
> 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]...
> [---=| TOFU protection by t-prot: 22 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
