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

author
- PiotrLi
See also
- https://introcs.cs.princeton.edu/java/97data/FFT.java.html
license
- 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)