#2 Functions: FFT, IFFT, IWT, IWT2,...

open
None
5
2003-10-23
2002-10-17
No

: Anyways, since i needed to implement those for
my own ATPG project, i decided to write those
functions for blitz since i will be doing math
using
blitz.

There could be bugs.

Anyways, i have implemented
1-d vectorized FFT/IFFT (tested)
Shuffle, BitReverse, FFT, IFFT

1-d vectorized FWT/IFWT (wavelet) (tested)
Uphi, DownHi, UpLo, DownLo, FWT_PO,
IWT_PO

2-d FWT/IFWT (wavelet) (tested, but not
completely)
FWT2_PO, IWT2_PO

1-d vectorized FWT( walsh transform) (tested)
FWT, IFWT

"vectorized" matrix inverse (tested)
Inv
matrix svd (tested)
SVD , conditional SVD
matrix pesudo inverse (tested)
PInv

matrix LU decomposition (tested)
LU, LUBSB, LUInverse

matrix Cholesky decompostion (untested)
LL, LLBSB

HighPrecision Arthimetic( uncompiled)

and have applied a TEMPORARY patch to Range
just to
get the range to be divisible by the stride.
[Not Included]

And added a dummy Array::ascending(int) const to
get
tensor multiplication to compile for version 0.6.
[Not Included]

thanks, hopefully u will like it.
-suresh

Discussion

  • Suresh kumar

    Suresh kumar - 2002-10-17

    functions

     
  • Julian Cummings

    Julian Cummings - 2003-10-23
    • assigned_to: nobody --> julianc
     
  • Shakes

    Shakes - 2006-06-30

    Logged In: YES
    user_id=1335435

    Re-Post:
    The following is the code for getting FFT's on Blitz Arrays
    using the mature and highly recognised FFTW C Library.

    template <typename type, int size>
    void
    FourierTransform<type,size>::FFTW(Array<complex<double>,2>
    &field, Array<complex<double>,2> &fftOfField)
    {
    complex<double> *ptrField, *ptrFFT;
    fftw_plan plan;

    ptrField = field.data();
    ptrFFT = fftOfField.data();
    plan =
    fftw_plan_dft_2d(field.rows(),field.cols(),reinterpret_cast<fftw_complex*>(ptrField),reinterpret_cast<fftw_complex*>(ptrFFT),FFTW_FORWARD,FFTW_ESTIMATE);

    fftw_execute(plan);

    fftw_destroy_plan(plan);
    }

    I have tested and its does work. Note that this works
    because as FFTW points out:
    "C++ has its own complex<T> template class, defined in the
    standard <complex> header file.
    Reportedly, the C++ standards committee has recently agreed
    to mandate that the storage
    format used for this type be binary-compatible with the C99
    type, i.e. an array T[2] with
    consecutive real [0] and imaginary [1] parts. (See report
    WG21/N1388.) Although not
    part of the official standard as of this writing, the
    proposal stated that: “This solution has
    been tested with all current major implementations of the
    standard library and shown to
    be working.” To the extent that this is true, if you have a
    variable complex<double> *x,
    you can pass it directly to FFTW via
    reinterpret_cast<fftw_complex*>(x)."

    Also data reordering has to be done if DC is wanted in the
    centre (http://www.fftw.org/faq/section3.html#centerorigin).
    I have done it in the following:
    void
    FourierTransform<type,size>::centerDC(Array<complex<double>,2>
    &field, Array<complex<double>,2> &result)
    {
    Array<complex<double>,2> tmp(field.rows(),field.cols());

    for(int k = 0; k < field.cols(); k ++)
    for(int j = 0; j < field.rows(); j ++)
    tmp((j+field.rows()/2)%field.rows(),(k+field.cols()/2)%field.cols())
    = field(j,k);

    result = tmp;
    }

    Code Excerpt is taken from the FourierTransform Class in the
    up and coming Quantum Mechanics and Discrete Geometry
    Toolkit based on Blitz++ and Qt - qC++
    (qcplusplus.sourceforge.net). It currently supports BEC
    simulations and the Fourier class can be found in the
    Subversion area.

    Keep up the good work guys, more work on the tiny stuff
    would be appreciated :)

    Cheers
    Shakes

     

Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks