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).

Turbo Prolog Index     Home Page