lambdaspeech :: lispology4

2 | list | login |load

| |

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

_h2 the Fast Fourier Transform {sup '{lambda talk} arrays}         
     
_img https://www.ritchievink.com/img/post-5-fft/fig_2.png

{center {h1 {code X(N) = Σ{sub n=0}{sup N-1} 1/n.sin(2Πn/N)}}}

_p We translate the JS code of the page [[JS_FFT|?view=lispology]] to build the Fast Fourier Transform algorithm on a few '{lambda talk} user defined functions.

_h3 1) the {code Cfft} function (direct & inverse)

{pre

'{def Cfft
 {lambda {:x :s}
  {if {> {#.length :x} 1}
   then {Ccombine :x :s 
         {Cfft {#.evens :x} :s} 
         {Cfft {#.odds  :x} :s} 0 {#.length :x}} 
   else :x
}}}
-> {def Cfft
 {lambda {:x :s}
  {if {> {#.length :x} 1}
   then {Ccombine :x :s
         {Cfft {#.evens :x} :s}
         {Cfft {#.odds  :x} :s} 0 {#.length :x}} 
   else :x
}}}

'{def Ccombine
 {lambda {:x :s :ev :od :k :N}
  {if {< :k {/ :N 2}}
   then {Ccombine {Cplusminus :x :s :ev :od :k :N}
                   :s :ev :od {+ :k 1} :N}
   else :x}}}
-> {def Ccombine
 {lambda {:x :s :ev :od :k :N}
  {if {< :k {/ :N 2}}
   then {Ccombine {Cplusminus :x :s :ev :od :k :N} :s :ev :od {+ :k 1} :N}
   else :x}}}

'{def Cplusminus
 {lambda {:x :s :ev :od :k :N}
  {let { {:x :x} {:k :k} {:N :N}
         {:evk {#.get :ev :k}}
         {:root {CUroot :s :od :k :N}} }
       {#.set! {#.set! :x {+ :k {/ :N 2}} {Csub :evk :root} } 
                                       :k {Cadd :evk :root} } }}} 
-> {def Cplusminus
 {lambda {:x :s :ev :od :k :N}
  {let { {:x :x} {:k :k} {:N :N}
         {:evk {#.get :ev :k}}
         {:root {CUroot :s :od :k :N}} }
       {#.set! {#.set! :x {+ :k {/ :N 2}} {Csub :evk :root} } 
                                       :k {Cadd :evk :root} } }}}

'{def CUroot                // x(k).exp(-/+i2πk/N)
 {lambda {:s :od :k :N}
  {Cmul {#.get :od :k} {Cexp {Cnew 0 {/ {* -2 :s {PI} :k} :N}}}} }}
-> {def CUroot
 {lambda {:s :od :k :N}
  {Cmul {#.get :od :k} {Cexp {Cnew 0 {/ {* -2 :s {PI} :k} :N}}}} }}

}

_h3 2) four functions on arrays

{pre
'{def #.evens
 {lambda {:x}
  {#.new {map {{lambda {:x :i} {#.get :x :i}} :x}
               {serie 0 {- {#.length :x} 1} 2}}} }}
-> {def #.evens
 {lambda {:x}
  {#.new {map {{lambda {:x :i} {#.get :x :i}} :x}
               {serie 0 {- {#.length :x} 1} 2}}} }}

'{def #.odds
 {lambda {:x}
  {#.new {map {{lambda {:x :i} {#.get :x :i}} :x}
               {serie 1 {#.length :x} 2}}} }}
-> {def #.odds
 {lambda {:x}
  {#.new {map {{lambda {:x :i} {#.get :x :i}} :x}
               {serie 1 {#.length :x} 2}}} }}

'{def #.sample2Complex
 {lambda {:x}
  {#.new {map {{lambda {:x :i} {Cnew {#.get :x :i} 0}} :x}
               {serie 0 {- {#.length :x} 1}}}} }}
-> {def #.sample2Complex
 {lambda {:x}
  {#.new {map {{lambda {:x :i} {Cnew {#.get :x :i} 0}} :x}
               {serie 0 {- {#.length :x} 1}}}} }}

'{def #.duplicate
 {lambda {:a}
  {#.slice :a 0 {#.length :a}} }}
-> {def #.duplicate
 {lambda {:a}
  {#.slice :a 0 {#.length :a}} }}

}

_h3 3) a small library for Cnumbers

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

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

_h3 4) testing

_p Following the example found in [[rosettacode.org/wiki/Fast_Fourier_transform|http://rosettacode.org/wiki/Fast_Fourier_transform#C.2B.2B]] we define a reduced sample of the square curve {code [1 1 1 1 0 0 0 0]}, translate it into Cnumbers, then apply the FFT {code '{Cfft x 1}} and the inverse FFT {code '{Cfft x -1}} to retrieve the initial sample.

{prewrap
'{def sample {#.new 1 1 1 1 0 0 0 0}}
-> {def sample {#.new 1 1 1 1 0 0 0 0}} 
'{sample}
-> {def sample {#.new 1 1 1 1 0 0 0 0}}

'{def x {#.sample2Complex {sample}}}
-> {def x {#.sample2Complex {sample}}} 
'{x}
-> {x}

'{def X {Cfft {#.duplicate {x}} 1}}
-> {def X {Cfft {#.duplicate {x}} 1}}
'{X}
-> {X}

'{def x' {Cfft {#.duplicate {X}} -1}}
-> {def x' {Cfft {#.duplicate {X}} -1}} 
'{x'}
-> {x'}

}

_p The last expression {code x'} is equivalent to the initial {code sample} where values are divided by 8 and values close to zero are displayed as 0.

_h3 5) finding the sines of the square curve

_p We sample the square curve with 128 levels, 64 at +1 and 64 at -1 and build an array of Cnumbers

{prewrap
'{def curve {#.new {map {lambda {_} 1} {serie 1 64}} 
                  {map {lambda {_} -1} {serie 1 64}} }}
-> {def curve {#.new
 {map {lambda {_} 1} {serie 1 64}} 
 {map {lambda {_} -1} {serie 1 64}} }}

'{curve}
-> {curve}
}

_p We apply the Fast Fourier Transform which returns an array of Cnumbers

{prewrap
'{def curve_fft {Cfft {#.sample2Complex {curve}} 1}}
-> {def curve_fft {Cfft {#.sample2Complex {curve}} 1}}

'{curve_fft {curve} 1}
-> {curve_fft {curve} 1}

}

_p We compute the norm of these Cnumbers and find the maximum value

{prewrap
'{def mods {map {{lambda {:x :i}
 {Cnorm {#.get :x :i}}} {curve_fft}} 
   {serie 0 {- {#.length {curve_fft}} 1}}} }
-> {def mods {map {{lambda {:x :i}
 {Cnorm {#.get :x :i}}} {curve_fft}} 
   {serie 0 {- {#.length {curve_fft}} 1}}} }

'{mods}
-> {mods}

'{def max_mods {max {mods}}}
-> {def max_mods {max {mods}}}
'{max_mods}
-> {max_mods}
}

_p We display the 20 first frequencies and amplitudes

{pre
'{map {lambda {:i} 
 {br}frequency: :i & amplitude: 
  {if {= {nth :i {mods}} 0}
   then . 
   else 1/{round {/ {max_mods} {nth :i {mods}}}}}}
    {serie 1 20}}
-> {map {lambda {:i} 
 {br}frequency: :i & amplitude: 
  {if {= {nth :i {mods}} 0}
   then . 
   else 1/{round {/ {max_mods} {nth :i {mods}}}}}}
    {serie 1 20}}
}

_p We have found the first sines of increasing odd frequency, [1,3,5,...,17,...], and decreasing amplitude, [1/1,1/3,1/3,...,1/17,...] leading to the square curve, with an exact precision from 1 to 17.

_p {i Alain Marty 2018/12/14}

{center {i Page computed in about 90ms on a recent laptop.}}

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