Cooley–Tukey FFT algorithm

   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(plcomplex).

 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    deal(X, Event, Odd),
   40
   41    ct_fft(Event,Q),
   42    ct_fft(Odd,R),
   43    Th is -2 * pi / N,
   44
   45    divide(Y, Y1, Y2),
   46
   47    % combine
   48    combine(Q,R,Th,0,Y1,Y2).
   49
   50combine([],[],_,_,_,_).
   51combine([Q|QT],[R|RT],Th,K,[Yk|Y1],[Ykn2|Y2]):-
   52    WkR iis R*exp(i*(K*Th)),
   53    Yk iis Q+WkR,
   54    Ykn2 iis Q-WkR,
   55    succ(K,KK),
   56    combine(QT,RT,Th,KK,Y1,Y2).
   57
   58deal([], [], []) :- !.
   59deal([A], [A], []) :- !.
   60deal([A,F,B,G,C,H,D,I,E,J|T],[A,B,C,D,E|L],[F,G,H,I,J|R]):-
   61    !,
   62    deal(T,L,R).
   63deal([A,F,B,G,C,H,D,I|T],[A,B,C,D|L],[F,G,H,I|R]):-
   64    !,
   65    deal(T,L,R).
   66deal([A,F,B,G,C,H|T],[A,B,C|L], [F,G,H|R]):-
   67    !,
   68    deal(T,L,R).
   69deal([A,F,B,G|T],[A,B|L], [F,G|R]):-
   70    !,
   71    deal(T,L,R).
   72deal([A,F|T],[A|L],[F|R]):-
   73    !,
   74    deal(T,L,R).
   75
   76divide(A,B,C):-
   77    divide_(A, B, A, C).
   78
   79divide_([], [], C, C).
   80divide_([_,_|T],[X|L],[X|LL], C):-
   81    divide_(T, L, LL, C).

 ct_ifft(+X:list, -Y:list) is det

Compute the inverse FFT of X, assuming its length is a power of 2.

   87ct_ifft([],[]):-!.
   88ct_ifft(X,Y):-
   89    length(X,N),
   90
   91    % take conjugate
   92    conjugate(X,Y1),
   93
   94    % compute forward FFT
   95    ct_fft(Y1,Y2),
   96
   97    % take conjugate again
   98    conjugate(Y2,Y3),
   99
  100    % divide by n
  101    maplist({N}/[AA,BB]>>(BB iis AA/N),Y3,Y),!.
  102
  103conjugate([],[]):-!.
  104conjugate([A|AT],[B|BT]):-
  105    B iis conjugate(A),
  106    conjugate(AT,BT),!.

 ct_cconvolve(+X:list, +Y:list, -Z:list) is det

Compute the circular convolution of X and Y.

  112ct_cconvolve(X,Y,Z):-
  113
  114    % compute FFT of each sequence
  115    ct_fft(X,A),
  116    ct_fft(Y,B),
  117
  118    % point-wise multiply
  119    maplist([AA,BB,CC]>>(CC iis AA*BB),A,B,C),
  120
  121    % compute inverse FFT
  122    ct_ifft(C,Z).

 ct_convolve(+X:list, +Y:list, -Z:list) is det

Compute the linear convolution of Y and Y.

  128ct_convolve(X,Y,Z):-
  129    length(X,N),
  130    N2 is N*2,
  131    length(A,N2),
  132    length(B,N2),
  133    length(ZZ,N),
  134    maplist([0]>>true,ZZ),
  135    append(X,ZZ,A),
  136    append(Y,ZZ,B),
  137    ct_cconvolve(A,B,Z)