annotate yetilab/transform/fft.yeti @ 271:c206de7c3018

Faster sparseProductLeft. Similar code would work for other sparseProducts

module yetilab.transform.fft;

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

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

load yetilab.vector.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 = new double[n*2-2];
    for [0..(n-1)*2-1] do i:
        ix = int (i/2);
        v[i] :=
            if i == ix*2 then
                complex.real arr[ix]
            else 
                complex.imaginary arr[ix] 
            fi;
    done;
    v[1] := complex.real arr[n-1];
    v);

//!!! doc: n is supplied separately from the input vector to support partial evaluation
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,
}