DFT | JS_FFT | LT_FFT ^{(arrays | lists)}

the Fast Fourier Transform ^{{lambda talk} lists}

1) the fft function (direct & inverse)

1.1) Following the standard FFT algorithm

{def fft   
 {lambda {:s :x}
  {if {= {list.length :x} 1}
   then :x
   else {let { {:s :s}
               {:ev {fft :s {evens :x}} }
               {:od {fft :s {odds  :x}} } }
        {let { {:ev :ev} {:t {rotate :s :od 0 {list.length :od}}} }
        {list.append {list.map Cadd :ev :t}
                     {list.map Csub :ev :t}} }}}}}

which can be rewritten without the let syntaxic sugar

{def fft     
 {lambda {:s :x} 
  {if {= {list.length :x} 1} 
   then :x
   else {{lambda {:s :ev :od}
          {{lambda {:ev :t}  
            {list.append {list.map Cadd :ev :t}
                         {list.map Csub :ev :t}} 
           } :ev {rotate :s :od 0 {list.length :od}}} 
         } :s {fft :s {evens :x}} {fft :s {odds :x}} }}}}
-> fft
  
{def rotate 
 {lambda {:s :f :k :N}
  {if {list.null? :f}
   then nil
   else {cons {Cmul {car :f} 
                    {Cexp {Cnew 0 {/ {* :s {PI} :k} :N}}}}
              {rotate :s {cdr :f} {+ :k 1} :N}}}}}
-> rotate

Note that the last version is the fastest one and could be rewritten without the {if then else} special form and reach the λ-calculus level with a small set of functions acting on words.

1.2) an alternate code

Following the clever and more Lisplike code written by Prasenjit Saha.

{def fft
 {lambda {:s :f}
  {if {= {list.length :f} 1}
   then :f
   else {combine :s {fft :s {evens :f}}
                    {fft :s {odds :f}}} }}}
{def combine
 {lambda {:s :ev :od}
  {plusminus :ev {rotate :s :od 0 {list.length :od}}}}}

{def plusminus
 {lambda {:a :b}
  {list.append {list.map Cadd :a :b}
               {list.map Csub :a :b}}}}

{def rotate 
 {lambda {:s :f :k :N}
  {if {list.null? :f}
   then nil
   else {cons {Cmul {car :f} 
                    {Cexp {Cnew 0 {/ {* :s {PI} :k} :N}}}}
              {rotate :s {cdr :f} {+ :k 1} :N}}}}}

2) adding a small set of list functions

We add to the existing {lambda talk}'s list primitives a small set of functions required by the function fft.

{def evens
 {lambda {:l}
  {if {list.null? :l}
   then nil
   else {cons {car :l} {evens {cdr {cdr :l}}}}}}}
-> evens

{def odds
 {lambda {:l}
  {if {list.null? {cdr :l}}
   then nil
   else {cons {car {cdr :l}} {odds {cdr {cdr :l}}}}}}}
-> odds

{def list.map
 {def list.map.r 
  {lambda {:f :a :b :c}
   {if {list.null? :a}
    then :c
    else {list.map.r :f {cdr :a} {cdr :b} 
                        {cons {:f {car :a} {car :b}} :c}} }}}
 {lambda {:f :a :b}
  {list.map.r :f {list.reverse :a} {list.reverse :b} nil}}}
-> list.map

3) adding a small set of complex functions

{lambda talk} has no primitive functions working on complex numbers. We add the minimal set required by the function fft.

{def Cnew
 {lambda {:x :y}
  {cons :x :y} }} 
-> Cnew

{def Cnorm
 {lambda {:c}
  {sqrt {+ {* {car :c} {car :c}}
           {* {cdr :c} {cdr :c}}}} }}
-> Cnorm

{def Cadd
 {lambda {:x :y}
  {cons {+ {car :x} {car :y}}
        {+ {cdr :x} {cdr :y}}} }}
-> Cadd

{def Csub
 {lambda {:x :y}
  {cons {- {car :x} {car :y}}
        {- {cdr :x} {cdr :y}}} }}
-> Csub

{def Cmul
 {lambda {:x :y}
  {cons {- {* {car :x} {car :y}} {* {cdr :x} {cdr :y}}}
        {+ {* {car :x} {cdr :y}} {* {cdr :x} {car :y}}}} }}
-> Cmul

{def Cexp
  {lambda {:x}
   {cons {* {exp {car :x}} {cos {cdr :x}}}
         {* {exp {car :x}} {sin {cdr :x}}}} }}
-> Cexp

{def Clist
 {lambda {:s}
  {list.new {map {lambda {:i} {cons :i 0}} :s}}}}
-> Clist

4) testing

Applying successively the direct fft function and its inverse on such a sample (1 -6 2 4) where numbers have been promoted as complex.

{list.disp {fft -1 {fft 1 {Clist 1 -6 2 4 }}}}

-> (4 0) (-24 1.0762083040283459e-16) (8 0) (16 -1.0762083040283459e-16)

returns Cnumbers whose real part divided by the length is the initial sample.

5) harmonics of a square curve

Applying the fft function on a sample (a1 a2 a3 ... an) will return a list of Cnumbers. We will compute the norm of these Cnumbers and display a table of harmonics [frequency,amplitude].

{def norms
 {lambda {:l}
  {if {list.null? :l}
   then nil
   else {cons {Cnorm {car :l}} {norms {cdr :l}}}}}}
-> norms
 
{def harmonics.disp
 {def harmonics.disp.r
  {lambda {:l :k :max :n}
   {if {> :k :n}
    then 
    else {br}:k: {if {= {car :l} 0} then . else 1/{round {/ :max {car :l}}}}      
         {harmonics.disp.r {cdr :l} {+ :k 1} :max :n}}}}
 {lambda {:l :N} {harmonics.disp.r :l 0 {max {list.disp :l}} :N}}}
-> harmonics.disp

We want to find the harmonics of the square curve whose sample is made of 32 levels +1 and 32 levels -1.

{def test {list.new 
  {map {lambda {_} {Cnew  1 0}} {serie 1 32}}
  {map {lambda {_} {Cnew -1 0}} {serie 1 32}} }} 
-> test    

{list.disp {test}} 
-> (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0) (-1 0)

We apply the function fft

{def X {fft -1 {test}}} 
-> X = (0 0) (2.0000000000000013 -40.71093524997437) (0 0) (2.000000000000002 -13.482904810829975) (0 0) (1.9999999999999996 -7.984447567540169) (0 0) (2.0000000000000018 -5.589625544980954) (0 0) (1.9999999999999991 -4.228644715097282) (0 0) (2.0000000000000018 -3.3367984111670146) (0 0) (1.9999999999999987 -2.6966878269734402) (0 0) (2.000000000000003 -2.206659951466955) (0 0) (1.9999999999999987 -1.8126943380382947) (0 0) (2.0000000000000004 -1.4833010925440706) (0 0) (2 -1.1987538673638474) (0 0) (2.000000000000001 -0.9459295517826396) (0 0) (1.9999999999999996 -0.7156114426290485) (0 0) (2.000000000000001 -0.5009739203826098) (0 0) (2 -0.29667197507669485) (0 0) (2.0000000000000044 -0.09825369953893315) (0 0) (1.9999999999999996 0.09825369953893315) (0 0) (1.9999999999999996 0.29667197507669485) (0 0) (2 0.5009739203826111) (0 0) (2 0.715611442629049) (0 0) (2.0000000000000004 0.9459295517826396) (0 0) (1.9999999999999998 1.1987538673638487) (0 0) (2 1.4833010925440697) (0 0) (2.000000000000003 1.8126943380382965) (0 0) (2.0000000000000004 2.20665995146695) (0 0) (1.9999999999999991 2.6966878269734402) (0 0) (1.9999999999999996 3.336798411167014) (0 0) (1.9999999999999984 4.228644715097283) (0 0) (2.0000000000000013 5.589625544980953) (0 0) (1.9999999999999956 7.984447567540169) (0 0) (2 13.482904810829975) (0 0) (1.999999999999982 40.71093524997437)

We compute norms

{def N {norms {X}}} 
-> N = 0 40.76003249419222 0 13.630433673874874 0 8.231124039813645 0 5.9366584652566665 0 4.6777597337315395 0 3.890272951447844 0 3.3573985816621073 0 2.9781450840091988 0 2.6992333658200223 0 2.4900164921426597 0 2.3317398728245355 0 2.212415584137779 0 2.1241703643591365 0 2.0617892396904987 0 2.0218838395901746 0 2.0024119929407895 0 2.0024119929407846 0 2.021883839590174 0 2.061789239690498 0 2.124170364359137 0 2.2124155841377786 0 2.331739872824536 0 2.490016492142659 0 2.6992333658200267 0 2.978145084009193 0 3.3573985816621077 0 3.8902729514478427 0 4.67775973373154 0 5.936658465256665 0 8.231124039813645 0 13.630433673874874 0 40.76003249419222

We display harmonics

{harmonics.disp {N} 15} 
->

0: . 

1: 1/1 

2: . 

3: 1/3 

4: . 

5: 1/5 

6: . 

7: 1/7 

8: . 

9: 1/9 

10: . 

11: 1/10 

12: . 

13: 1/12 

14: . 

15: 1/14

We have found the first odd sines [1,1/1] [3,1/3] [5,1/5] [7,1/7] [9,1/9].

Increasing the size of sample would give some of the following sines:

And so on.

Next step: reach the λ-calculus level.

Your feedback is welcome in the wiki's agora or in HN. See also rosetta.

Alain Marty 2018/12/(21 edit 27)

Edit: Have a look at this amazing work http://www.jezzamon.com/fourier.