Chris@41: 
Chris@93: module yetilab.transform.fft;
Chris@41: 
Chris@41: import edu.emory.mathcs.jtransforms.fft: DoubleFFT_1D;
Chris@41: 
Chris@222: vec = load yetilab.vector.vector;
Chris@222: complex = load yetilab.vector.complex;
Chris@41: 
Chris@222: load yetilab.vector.complextype;
Chris@213: 
Chris@213: packedToComplex len p is number -> ~double[] -> array<cplx> =
Chris@213:    (n = len / 2;
Chris@46:     array
Chris@46:        (map do i:
Chris@46:             re = if i == n then p[1] else p[i*2] fi;
Chris@46:             im = if i == 0 or i == n then 0 else p[i*2+1] fi;
Chris@46:             complex.complex re im;
Chris@46:         done [0..n]));
Chris@41: 
Chris@48: complexToPacked arr =
Chris@48:    (n = length arr;
Chris@221:     v = new double[n*2-2];
Chris@221:     for [0..(n-1)*2-1] do i:
Chris@221:         ix = int (i/2);
Chris@221:         v[i] :=
Chris@48:             if i == ix*2 then
Chris@48:                 complex.real arr[ix]
Chris@48:             else 
Chris@48:                 complex.imaginary arr[ix] 
Chris@48:             fi;
Chris@221:     done;
Chris@48:     v[1] := complex.real arr[n-1];
Chris@48:     v);
Chris@48: 
Chris@48: realForward n = 
Chris@41:    (d = new DoubleFFT_1D(n);
Chris@41:     do bl:
Chris@213:         v = vec.primitive bl;
Chris@41:         d#realForward(v);
Chris@213:         packedToComplex (vec.length bl) v;
Chris@48:     done);
Chris@48: 
Chris@48: realInverse n = 
Chris@48:    (d = new DoubleFFT_1D(n);
Chris@48:     do arr:
Chris@48:         v = complexToPacked arr;
Chris@48:         d#realInverse(v, true);
Chris@213:         vec.vector v;
Chris@41:     done);
Chris@41: 
Chris@41: {
Chris@48: realForward,
Chris@48: realInverse,
Chris@41: }
Chris@41: