lambdaspeech

2 | list | login |load

| |

_h3 DFT | [[JS_FFT|?view=lispology]] | [[LT_FFT|?view=lispology4]]

_h2 the Discrete Fourier Transform

_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)}}}

_p We progressively approximate a {b square curve} as the sum of sines of increasing odd frequencies and decreasing amplitudes: [ [1,1/1] [3,1/3] [5,1/5] [7,1/7] [9,1/9] ... ].

{pre
'{def CURVE 
 {lambda {:k :t}
  {- :t 180}
  {* 100 {+ {map {{lambda {:t :k}
                   {* {/ 1 :k} {sin {* {/ {PI} 180} :k :t}}} } :t} 
                 {serie 1 :k 2}}} }}}
-> {def CURVE
 {lambda {:k :t}
  {- :t 180}
  {* 140 {+ {map {{lambda {:t :k}
                   {* {/ 1 :k} {sin {* {/ {PI} 180} :k :t}}} }
                  :t} {serie 1 :k 2}}} }}}
}

_p We plot the first five curves

{pre
'{svg {@ width="600px" height="300px" style="background:#eee"}
 {g {AXES 600 300}
  {polyline {@ points="-180 0 -180 110 0 110 0 -110 180 -110 180 0" {stroke #888 5"}}}
  {polyline {@ points="{map {CURVE 1} {serie 0 360 2}}" {stroke #222 1}}}
  {polyline {@ points="{map {CURVE 3} {serie 0 360 2}}" {stroke #444 1}}}
  {polyline {@ points="{map {CURVE 5} {serie 0 360 2}}" {stroke #444 1}}}
  {polyline {@ points="{map {CURVE 7} {serie 0 360 2}}" {stroke #444 1}}}
  {polyline {@ points="{map {CURVE 9} {serie 0 360 2}}" {stroke #f00 3}}} 
}}
}

{svg {@ width="600px" height="300px" style="background:#eee"}
 {g {AXES 600 300}
  {polyline {@ points="-180 0 -180 110 0 110 0 -110 180 -110 180 0" {stroke #888 5"}}}
  {polyline {@ points="{map {CURVE 2} {serie 0 360 2}}" {stroke #222 1}}}
  {polyline {@ points="{map {CURVE 3} {serie 0 360 2}}" {stroke #444 1}}}
  {polyline {@ points="{map {CURVE 5} {serie 0 360 2}}" {stroke #444 1}}}
  {polyline {@ points="{map {CURVE 7} {serie 0 360 2}}" {stroke #444 1}}}
  {polyline {@ points="{map {CURVE 9} {serie 0 360 2}}" {stroke #f00 3}}} 
}}

_p The square curve is the limit of this process. Such a curve, two horizontal segments at levels +100 and -100, can be sampled as, say 32 levels at +100 and 32 levels at -100

{prewrap
'{def sample 
 {list.new {map {lambda {_}  100} {serie 1 32}} 
           {map {lambda {_} -100} {serie 1 32}} }}
-> {def sample {list.new {map {lambda {_} 100} {serie 1 32}} {map {lambda {_} -100} {serie 1 32}} }} 
'{list.disp {sample}} -> 
{list.disp {sample}}
}

_p Given this sample we want to find back the frequencies and amplitudes. Following [[http://practicalcryptography.com|http://practicalcryptography.com/miscellaneous/machine-learning/intuitive-guide-discrete-fourier-transform/]] or/and [[https://www.nayuki.io|https://www.nayuki.io/page/how-to-implement-the-discrete-fourier-transform]], we will use the {b Discrete Fourier Transform} (DFT) which is an algorithm that takes a sample of some curve and determines its {b frequency spectrum}. The DFT algorithm computes the vectorial product of a sequence of N values and the roots of unity given {code e{sup i2π*n}}. Its expression is usualy given like this

{center {pre
X(x,k) = {SIGMA}x(n)*e{sup -i2π*kn/N}
}}

_p but is better rewritten as a double sum of trigonometric functions

{center {pre
X(x,k) = {SIGMA}x(n)*cos(2π*kn/N)-i{SIGMA}x(n)*sin(2π*kn/N)
}}

_p In fact we will build a {b dft(x,k)} function taking as input the sample {b x(n)} and a frequency {b k} and returning the norm of {b X(k)}, {b √(X.real{sup 2}+X.imag{sup 2})}.

{pre
'{def dft
 {def dft.r {lambda {:x :k :X :Y :n}
   {if {= :n 0}
    then {round {sqrt {+ {* :X :X} {* :Y :Y}}}}
    else {dft.r {cdr :x} 
                :k
                {+ :X {* {car :x} {cos {* :k :n}}}}   
                {+ :Y {* {car :x} {sin {* :k :n}}}} 
                {- :n 1}} 
 }}}
 {lambda {:x :k}
  {dft.r :x 
         {/ {* 2 {PI} :k} {list.length :x}} 
         0 
         0 
         {- {list.length :x} 1}} 
}} 
-> {def dft
 {def dft.r {lambda {:x :k :X :Y :n}
   {if {= :n 0}
    then {round {sqrt {+ {* :X :X} {* :Y :Y}}}}
    else {dft.r {cdr :x} 
                :k
                {+ :X {* {car :x} {cos {* :k :n}}}}   
                {+ :Y {* {car :x} {sin {* :k :n}}}} 
                {- :n 1}} 
 }}}
 {lambda {:x :k}
  {dft.r :x 
         {/ {* 2 {PI} :k} {list.length :x}} 
         0 
         0 
         {- {list.length :x} 1}} 
}}
}

_p We compute correlations for frequencies in range [1,15]

{pre 
'{def dfts {map {dft {sample}} {serie 1 15}}} 
-> {def dfts {map {dft {sample}} {serie 1 15}}}

'{dfts}
-> {dfts}
}

_p We find the maximum dft value and display frequencies and amplitudes

{pre
'{def dft.max {max {dfts}}} 
-> {def dft.max {max {dfts}}} = {dft.max}

'{map {lambda {:k}
       {br}{+ :k 1}:
           1/{round {/ {dft.max} {nth :k {dfts}}}}}
     {serie 0 14}}
->
}

{center {pre {@ style="text-align:left; width:100px;"}
{map {lambda {:k} {br}{+ :k 1}: 1/{round {/ {dft.max} {nth :k {dfts}}}}}
     {serie 0 14}}
}}

_p Or in other words

{pre
'{svg {@ width="600px" height="600px" style="background:#eee"}
 {g {AXES 600 600}
  {polyline {@ points="{map {lambda {:k}
                               {* 20 {+ :k 1}}
                               {* 200 {/ {nth :k {dfts}} {dft.max} }}} 
                            {serie 0 14}}"
            {stroke #00f 1}}} 
}}
}

{svg {@ width="600px" height="600px" style="background:#eee"}
 {g {AXES 600 600}
  {polyline {@ points="{map {lambda {:k} {* 20 {+ :k 1}} {* 200 {/ {nth :k {dfts}} {dft.max} }}} {serie 0 14}}" {stroke #00f 1}}}
}}

_p We have got back the five first increasing odd frequencies and decreasing amplitudes and some others leading to the square curve.

{div {@ style="font:italic 1.5em courier; text-align:center;"} C(t) = Σ{sub n=0}{sup ∞} 1/n*sin(π/180*n*t) }

_p {i Alain Marty 2018/12/07}

{{hide}
{def SIGMA
 {table {@ style="display:inline-block; vertical-align:middle;"}
  {tr {td N-1}}
  {tr {td {@ style="font-size:3.0em; line-height:0.5em;"}Σ}}
  {tr {td n=0}}
}}

{def AXES 
 {lambda {:w :h}
  {@ transform="translate({/ :w 2},{/ :h 2}) scale(1,-1)"} 
  {line {@ x1="-{/ :w 2}:w" y1="0" 
           x2="{/ :w 2}" y2="0"
           stroke="red" fill="transparent"}}
  {line {@ x1="0" y1="-{/ :h 2}" 
           x2="0" y2="{/ :h 2}"
           stroke="green" fill="transparent"}} 
}}
 
{def stroke
 {lambda {:col :w}
 stroke=":col" fill="transparent" stroke-width=":w" }}
}

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