#### Fast Fourier Transform

check_determ
domains
cplx  = c(real,real)
cplxs = cplx*
int   = integer
ints  = int*
reals = real*
predicates
fft(int,cplxs,cplxs)
fft_1(cplxs,cplxs,cplx,cplx,cplxs,cplxs)
evens_and_odds(reals,reals,reals)
evens_and_odds(ints,ints,ints)
evens_and_odds(cplxs,cplxs,cplxs)
sum(cplx,cplx,cplx)
product(cplx,cplx,cplx)
w(int,cplx)
fff(int,cplxs,cplxs)
left_and_right(int,cplxs,cplxs,cplxs)
fff_1(cplxs,cplxs,cplx,cplx,cplxs,cplxs)
diff(cplx,cplx,cplx)
inverse_fft(int,cplxs,cplxs)
inverse_fff(int,cplxs,cplxs)
conjugates(cplxs,cplxs)
scale(cplxs,int,cplxs)
clauses

/*
* Decimation-in-Time (Cooley-Tukey) Algorithm
*/

/* fft(N, F, Ft) is true if the list Ft contains the Fourier transform of  */
/*   the N -- a power of two -- samples in the list F.  Each sample is a   */
/*   complex number represented by c(RealPart, ImaginaryPart).             */
fft(1, F, F):-
!.
fft(N, F, Ft):-
N > 1,
N1 = N div 2,
evens_and_odds(F, E, O),
fft(N1, E, Et),
fft(N1, O, Ot),
w(1, W1),
w(2, W2),
w(N, Wn),
fft_1(Et, Ot, W2, Wn, Gt, []),
fft_1(Et, Ot, W1, Wn, Ft, Gt).

fft_1([], [], _, _, Ft, Ft):-!.
fft_1([E|Et], [O|Ot], Wk, Wn, [F|Ft], Fu):-
product(O, Wk, Temp),
sum(E, Temp, F),
product(Wk, Wn, Wk1),
fft_1(Et, Ot, Wk1, Wn, Ft, Fu).

/* evens_and_odds(Xs, Evens, Odds) is true if Evens is the list of the     */
/*   even-positioned elements of the list Xs, and Odds is the list of the  */
/*   odd-positioned elements of the list Xs, where Xs has an even number   */
/*   of elements, and the first element of Xs is considered to be at an    */
/*   even position.                                                        */
evens_and_odds([], [], []).
evens_and_odds([Even,Odd|Xs], [Even|Evens], [Odd|Odds]):-
evens_and_odds(Xs, Evens, Odds).

/* sum(A, B, C) is true if C is the sum of the complex numbers A and B.    */
sum(c(Ra,Ia), c(Rb,Ib), c(Rc,Ic)):-
Rc = Ra + Rb,
Ic = Ia + Ib.

/* product(A, B, C) is true if C is the product of the complex numbers A   */
/*   and B.                                                                */
product(c(Ra,Ia), c(Rb,Ib), c(Rc,Ic)):-
Rc = Ra*Rb - Ia*Ib,
Ic = Ra*Ib + Ia*Rb.

/* w(N, c(Cosine, Sine)) is true if Cosine = cos(2 * pi / N) and           */
/*   Sine = sin(2 * pi / N), where N is a power of two between 1 and 16384 */
/*   inclusive.                                                            */
w(    1, c( 1.0, 0.0)).
w(    2, c(-1.0, 0.0)).
w(    4, c( 0.0, 1.0)).
w(    8, c( 0.70710678119, 0.70710678119)).
w(   16, c( 0.92387953251, 0.38268343237)).
w(   32, c( 0.9807852804,  0.19509032202)).
w(   64, c( 0.99518472667, 0.09801714033)).
w(  128, c( 0.99879545621, 0.049067674327)).
w(  256, c( 0.9996988187,  0.024541228523)).
w(  512, c( 0.99992470184, 0.012271538286)).
w( 1024, c( 0.99998117528, 0.0061358846492)).
w( 2048, c( 0.99999529381, 0.003067956763)).
w( 4096, c( 0.99999882345, 0.0015339801863)).
w( 8192, c( 0.99999970586, 0.00076699031874)).
w(16384, c( 0.99999992647, 0.00038349516882)).

/*
* Decimation-in-Frequency (Sande-Tukey) Algorithm
*/

/* fff(N, F, Ft) is true if the list Ft contains the Fourier transform of  */
/*   the N -- a power of two -- samples in the list F.  Each sample is a   */
/*   complex number represented by c(RealPart, ImaginaryPart).             */
fff(1, F, F):-
!.
fff(N, F, Ft):-
N > 1,
N1 = N div 2,
left_and_right(N1, F, Left, Right),
w(1, W1),
w(N, Wn),
fff_1(Left, Right, W1, Wn, Sums, Diffs),
fff(N1, Sums, Et),
fff(N1, Diffs, Ot),
evens_and_odds(Ft, Et, Ot).

fff_1([], [], _, _, [], []).
fff_1([A|As], [B|Bs], Wn, W, [C|Cs], [D|Ds]):-
sum(A, B, C),
diff(A, B, D0),
product(D0, Wn, D),
product(Wn, W, Wn1),
fff_1(As, Bs, Wn1, W, Cs, Ds).

/* left_and_right(N, As, Bs, Cs) is true if Bs contains the first N        */
/*   elements of As, and Cs contains the remaining elements of As.         */
left_and_right(0, As, [], As):-!.
left_and_right(N, [A|As], [A|Bs], Cs):-
N1 = N - 1,
left_and_right(N1, As, Bs, Cs).

/* diff(A, B, C) is true if C is the difference between the complex        */
/*   numbers A and B.                                                      */
diff(c(Ra,Ia), c(Rb,Ib), c(Rc,Ic)):-
Rc = Ra - Rb,
Ic = Ia - Ib.

/*
* Inverse Fast Fourier Transforms
*/

/* inverse_fft(N, Ft, F) is true if the list F contains the inverse        */
/*   Fourier transform of the N -- a power of two -- samples in the list   */
/*   Ft.  Each sample is a complex number represented by c(Real, Imag).    */
inverse_fft(N, Ft, F):-
conjugates(Ft, Fc),
fft(N, Fc, Fct),
conjugates(Fct, Fn),
scale(Fn, N, F).

/* inverse_fff(N, Ft, F) is true if the list F contains the inverse        */
/*   Fourier transform of the N -- a power of two -- samples in the list   */
/*   Ft.  Each sample is a complex number represented by c(Real, Imag).    */
inverse_fff(N, Ft, F):-
conjugates(Ft, Fc),
fff(N, Fc, Fct),
conjugates(Fct, Fn),
scale(Fn, N, F).

/* conjugates(As, Bs) is true if the elements in the list Bs are the       */
/*   conjugates of the elements -- complex numbers -- in the list As.      */
conjugates([], []).
conjugates([c(Ra,Ia)|As], [c(Ra,Ib)|Bs]):-
Ib = -Ia,
conjugates(As, Bs).

/* scale(As, N, Bs) is true if the elements in the list Bs are the         */
/*   elements -- complex numbers -- in the list As divided by N.           */
scale([], _, []).
scale([c(Ra,Ia)|As], N, [c(Rb,Ib)|Bs]):-
Rb = Ra / N,
Ib = Ia / N,
scale(As, N, Bs).