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