# 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 -----------

• assumes length is a power of 2
• not the most memory efficient algorithm (because it uses an object type for representing complex numbers and because it re-allocates memory for the subarray, instead of doing in-place or reusing a single temporary array)
author
- PiotrLi
- https://introcs.cs.princeton.edu/java/97data/FFT.java.html
- GPL /
```   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).```
ct_fft(+X:list, -Y:list) is det
Compute the FFT of X, assuming its length is a power of 2.
```   33ct_fft([],[]):-!.
34ct_fft([X],[X]):-!.
35ct_fft(X,Y):-
36    length(X,N),
37    length(Y,N),
38
39    % fft of even terms
40    every_second(X,Event),
41
42    % fft of odd terms
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    % combine
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),!.```
ct_ifft(+X:list, -Y:list) is det
Compute the inverse FFT of X, assuming its length is a power of 2.
```   82ct_ifft([],[]):-!.
83ct_ifft(X,Y):-
84    length(X,N),
85
86    % take conjugate
87    conjugate(X,Y1),
88
89    % compute forward FFT
90    ct_fft(Y1,Y2),
91
92    % take conjugate again
93    conjugate(Y2,Y3),
94
95    % divide by n
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),!.```
ct_cconvolve(+X:list, +Y:list, -Z:list) is det
Compute the circular convolution of X and Y.
```  107ct_cconvolve(X,Y,Z):-
108
109    % compute FFT of each sequence
110    ct_fft(X,A),
111    ct_fft(Y,B),
112
113    % point-wise multiply
114    maplist([AA,BB,CC]>>(CC iis AA*BB),A,B,C),
115
116    % compute inverse FFT
117    ct_ifft(C,Z).```
ct_convolve(+X:list, +Y:list, -Z:list) is det
Compute the linear convolution of Y and Y.
```  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)```