view yetilab/transform/fft.yeti @ 213:d841bfc143dd matrix_opaque_immutable

Update complex, fft

module yetilab.transform.fft;

import edu.emory.mathcs.jtransforms.fft: DoubleFFT_1D;

vec = load yetilab.block.vector;
complex = load yetilab.block.complex;

load yetilab.block.complextype;

packedToComplex len p is number -> ~double[] -> array<cplx> =
   (n = len / 2;
    array
       (map do i:
            re = if i == n then p[1] else p[i*2] fi;
            im = if i == 0 or i == n then 0 else p[i*2+1] fi;
            complex.complex re im;
        done [0..n]));

complexToPacked arr =
   (n = length arr;
    v = (map do i:
            ix = int (i/2);
            if i == ix*2 then
                complex.real arr[ix]
            else 
                complex.imaginary arr[ix] 
            fi;
            done [0..(n-1)*2-1]) as ~double[];
    v[1] := complex.real arr[n-1];
    v);

realForward n = 
   (d = new DoubleFFT_1D(n);
    do bl:
        v = vec.primitive bl;
        d#realForward(v);
        packedToComplex (vec.length bl) v;
    done);

realInverse n = 
   (d = new DoubleFFT_1D(n);
    do arr:
        v = complexToPacked arr;
        d#realInverse(v, true);
        vec.vector v;
    done);

{
realForward,
realInverse,
}