lambdaspeech

lambdaspeech :: zorg

1 | list | login |load

| |

_h3 [[DFT]] | [[JS_FFT|?view=lispology]] | LT_FFT {sup ([[arrays|?view=lispology4]] | lists)}

_h2 the Fast Fourier Transform {sup '{lambda talk} lists}

_img https://www.ritchievink.com/img/post-5-fft/fig_2.png

{center {h1 {code X(N) = Σ{sub n=0}{sup {@ style="margin-left:-40px;"} N-1}  1/n.sin(2Πn/N)}}}

_h2 1) the {b fft} function (direct & inverse)

_h3 1.1) Following the standard FFT algorithm

{pre    
'{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}} }}}}} 
}

_p which can be rewritten without the let syntaxic sugar

{pre
'{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}} }}}}
-> {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}} }}}}
  
'{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}}}}}
-> {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}}}}}
}

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

_h3 1.2) an alternate code

_p Following the clever and more Lisplike code written by [[Prasenjit Saha|https://www.physik.uzh.ch/~psaha/mus/fourlisp.php]].

{pre
'{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}}}}}
}

_h2 2) adding a small set of list functions

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

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

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

'{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}}}
-> {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}}}
}

_h2 3) adding a small set of complex functions

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

{pre
;; {def isC {lambda {:c} {pair? :c}}}
'{def Cnew
 {lambda {:x :y}
  {cons :x :y} }} 
-> {def Cnew {lambda {:x :y} {cons :x :y} }}

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

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

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

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

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

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

}

_h2 4) testing

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

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

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

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

_h2 5) harmonics of a square curve

_p Applying the {code fft} function on a sample {code (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].

{prewrap
'{def norms
 {lambda {:l}
  {if {list.null? :l}
   then nil
   else {cons {Cnorm {car :l}} {norms {cdr :l}}}}}}
-> {def norms
 {lambda {:l}
  {if {list.null? :l}
   then nil
   else {cons {Cnorm {car :l}} {norms {cdr :l}}}}}}
 
'{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}}}
-> {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}}}

}

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

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

'{list.disp {test}} 
-> {list.disp {test}}
}

_p We apply the function {code fft}

{prewrap
'{def X {fft -1 {test}}} 
-> {def X {fft -1 {test}}} = {list.disp {X}} 
}

_p We compute norms

{prewrap
'{def N {norms {X}}} 
-> {def N {norms {X}}} = {list.disp {N}}
}

_p We display harmonics

{pre
'{harmonics.disp {N} 15} 
-> 
}
{center {pre {@ style="width:100px; text-align:left; background:#eee; padding:20px; box-shadow:0 0 8px #000;"}
{harmonics.disp {N} 15}
}}

_p We have found the first odd sines [1,1/1] [3,1/3] [5,1/5] [7,1/7]  [9,1/9]. 
_p Increasing the size of sample would give some of the following sines:

{center
{table {@ style="width:50%; text-align:center; background:#eee; box-shadow:0 0 8px #000;"}
{tr {td {b sample}} {td {b last sine}} {td {b time ms}} {td {b t{sub i+1}/t{sub i}}}}

{tr {td 8} {td [1,1/1]} {td 12} {td -}}
{tr {td 16} {td [3,1/3]} {td 17} {td {format {/ 17 12}}}}
{tr {td 32} {td [5,1/5]} {td 25} {td {format {/ 25 17}}}}
{tr {td 64} {td [9,1/9]} {td 45} {td {format {/ 45 25}}}}
{tr {td 128} {td [17,1/17]} {td 100} {td {format {/ 100 45}}}}
{tr {td 512} {td [43,1/43]} {td 658} {td {format {/ 658 100}}}}
{tr {td 1024} {td [67,1/67]} {td 1400} {td {format {/ 1400 658}}}}
{tr {td 2048} {td [107,1/107]} {td 4100} {td {format {/ 4100 1400}}}}
}}

_p And so on.

_p {b Next step}: reach the [[λ-calculus|?view=PLR]] level.

_p Your feedback is welcome in the wiki's [[agora]] or in [[HN|https://news.ycombinator.com/item?id=18761359]]. See also [[rosetta|http://www.rosettacode.org/wiki/Fast_Fourier_transform#lambdatalk]].

_p {i Alain Marty 2018/12/(21 edit 27)}

{blockquote {b Edit}: {i Have a look at this amazing work [[http://www.jezzamon.com/fourier|http://www.jezzamon.com/fourier/index.html]].}}

{{hide}
 {def format {lambda {:x} {/ {round {* 1000 :x}} 1000}}} 
}

{style ;;
@import url(https://fonts.googleapis.com/css?family=Quicksand);

#page_frame { width:620px; }

#page_content {
 font-family: Quicksand; font-size:1.0em;
 background:#ffe;
} 
}