20:-module(ct_fft,[
21 ct_fft/2,
22 ct_ifft/2,
23 ct_cconvolve/3,
24 ct_convolve/3
25 ]). 26
27:-use_module(complex).
33ct_fft([],[]):-!.
34ct_fft([X],[X]):-!.
35ct_fft(X,Y):-
36 length(X,N),
37 length(Y,N),
38
39 40 every_second(X,Event),
41
42 43 every_second([k|X],[k|Odd]),
44 ct_fft(Event,Q),
45 ct_fft(Odd,R),
46 Th is -2 * pi / N,
47 N2 is N/2,
48
49 50 combine(Q,R,Th,0,N2,Y),!.
51
52combine([],[],_,_,_,_):-!.
53combine([Q|QT],[R|RT],Th,K,N2,Y):-
54 Kth is K*Th,
55 Wk iis 1*exp(i*Kth),
56 WkR iis R*Wk,
57 Yk iis Q+WkR,
58 Ykn2 iis Q-WkR,
59 Kn2 is K + N2,
60 nth0(K,Y,Yk),
61 nth0(Kn2,Y,Ykn2),
62 succ(K,KK),
63 combine(QT,RT,Th,KK,N2,Y),!.
64
65every_second([],[]):-!.
66every_second([A],[A]):-!.
67every_second([A,_,B,_,C,_,D,_,E,_|T],[A,B,C,D,E|TT]):-
68 every_second(T,TT),!.
69every_second([A,_,B,_,C,_,D,_|T],[A,B,C,D|TT]):-
70 every_second(T,TT),!.
71every_second([A,_,B,_,C,_|T],[A,B,C|TT]):-
72 every_second(T,TT),!.
73every_second([A,_,B,_|T],[A,B|TT]):-
74 every_second(T,TT),!.
75every_second([A,_|T],[A|TT]):-
76 every_second(T,TT),!.
82ct_ifft([],[]):-!.
83ct_ifft(X,Y):-
84 length(X,N),
85
86 87 conjugate(X,Y1),
88
89 90 ct_fft(Y1,Y2),
91
92 93 conjugate(Y2,Y3),
94
95 96 maplist({N}/[AA,BB]>>(BB iis AA/N),Y3,Y),!.
97
98conjugate([],[]):-!.
99conjugate([A|AT],[B|BT]):-
100 B iis conjugate(A),
101 conjugate(AT,BT),!.
107ct_cconvolve(X,Y,Z):-
108
109 110 ct_fft(X,A),
111 ct_fft(Y,B),
112
113 114 maplist([AA,BB,CC]>>(CC iis AA*BB),A,B,C),
115
116 117 ct_ifft(C,Z).
123ct_convolve(X,Y,Z):-
124 length(X,N),
125 N2 is N*2,
126 length(A,N2),
127 length(B,N2),
128 length(ZZ,N),
129 maplist([0]>>true,ZZ),
130 append(X,ZZ,A),
131 append(Y,ZZ,B),
132 ct_cconvolve(A,B,Z)
Cooley–Tukey FFT algorithm
Compute the FFT and inverse FFT of a length n complex sequence using the radix 2 Cooley-Tukey algorithm.
Bare bones implementation that runs in O(n log n) time.
Limitations -----------