Fast Fourier Transform

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 is N // 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 the first element of Xs */
/*   is considered to be at an even position.                              */
evens_and_odds([], [], []).
evens_and_odds([X|Xs], [X|Ys], Zs):-
  evens_and_odds(Xs, Zs, Ys).

/* 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 is Ra + Rb,
  Ic is 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 is Ra*Rb - Ia*Ib,
  Ic is 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 is N // 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 is 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 is Ra - Rb,
  Ic is 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 is -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 is Ra / N,
  Ib is Ia / N,
  scale(As, N, Bs).

LPA Index     Home Page