lambdaway :: fft

8 | list | login |load

| |

_img data/FFT.png

_h1 FFT | [[fft]] | [[dft]] | [[zorg]]

{require lib_complex}

_ul [[../lambdaspeech/?view=FFT|../lambdaspeech/?view=FFT]]
_ul [[http://jeffe.cs.illinois.edu/teaching/algorithms/notes/A-fft.pdf|http://jeffe.cs.illinois.edu/teaching/algorithms/notes/A-fft.pdf]].

_h2 1) code

{pre
'{def A.fft   
 {lambda {:s :x}
  {if {= {A.length :x} 1}
   then :x
   else {let { {:s :s}
               {:ev {A.fft :s {A.evens :x}} }
               {:od {A.fft :s {A.odds  :x}} } }
               {let { {:ev :ev}
                      {:t {A.rotate :s :od 0 {A.length :od}}} }
                    {A.concat {A.apply C.add :ev :t}
                              {A.apply C.sub :ev :t}} }}}}}
-> {def A.fft   
 {lambda {:s :x}
  {if {= {A.length :x} 1}
   then :x
   else {let { {:s :s}
               {:ev {A.fft :s {A.evens :x}} }
               {:od {A.fft :s {A.odds  :x}} } }
        {let { {:ev :ev} {:t {A.rotate :s :od 0 {A.length :od}}} }
             {A.concat {A.apply C.add :ev :t}
                       {A.apply C.sub :ev :t}} }}}}}

'{def A.evens
 {def A.evens.r
  {lambda {:a :b}
   {if {A.empty? :a}
    then :b
    else {A.evens.r {A.rest {A.rest :a}}
                    {A.addlast! {A.first :a} :b}}}}}
 {lambda {:a}
  {A.evens.r :a {A.new}}}}
-> {def A.evens
 {def A.evens.r
  {lambda {:a :b}
   {if {A.empty? :a}
    then :b
    else {A.evens.r {A.rest {A.rest :a}}
                    {A.addlast! {A.first :a} :b}}}}}
 {lambda {:a}
  {A.evens.r :a {A.new}}}}

'{def A.odds
 {def A.odds.r
  {lambda {:a :b}
   {if {A.empty? {A.rest :a}} 
    then :b
    else {A.odds.r {A.rest {A.rest :a}}
                   {A.addlast! {A.first {A.rest :a}} :b}}}}}
 {lambda {:a}
  {A.odds.r :a {A.new}}}}
-> {def A.odds
 {def A.odds.r
  {lambda {:a :b}
   {if {A.empty? {A.rest :a}} 
    then :b
    else {A.odds.r {A.rest {A.rest :a}}
                   {A.addlast! {A.first {A.rest :a}} :b}}}}}
 {lambda {:a}
  {A.odds.r :a {A.new}}}}

'{def A.rotate
 {def A.rotate.r 
  {lambda {:s :f :k :N :b}
   {if {A.empty? :f} 
    then :b
    else {A.rotate.r :s {A.rest :f} {+ :k 1} :N
                     {A.addlast! {C.mul {A.first :f}
                                        {C.exp {C.new 0 {/ {* :s {PI} :k} :N}}}}
                                 :b}}}}}
 {lambda {:s :f :k :N} 
   {A.rotate.r :s :f :k :N {A.new}}}}
-> {def A.rotate
 {def A.rotate.r 
  {lambda {:s :f :k :N :b}
   {if {A.empty? :f} 
    then :b
    else {A.rotate.r :s {A.rest :f} {+ :k 1} :N
                     {A.addlast! {C.mul {A.first :f}
                                 {C.exp {C.new 0 {/ {* :s {PI} :k} :N}}}} :b}}}}}
 {lambda {:s :f :k :N} 
   {A.rotate.r :s :f :k :N {A.new}}}}

'{def A.apply
 {def A.apply.r 
  {lambda {:f :a :b :c}
   {if {A.empty? :a}
    then :c
    else {A.apply.r :f {A.rest :a} {A.rest :b} 
                       {A.addlast! {:f {A.first :a} {A.first :b}} :c}} }}}
 {lambda {:f :a :b}
  {A.apply.r :f :a :b {A.new}}}}
-> {def A.apply
 {def A.apply.r 
  {lambda {:f :a :b :c}
   {if {A.empty? :a}
    then :c
    else {A.apply.r :f {A.rest :a} {A.rest :b} 
                       {A.addlast! {:f {A.first :a} {A.first :b}} :c}} }}}
 {lambda {:f :a :b}
  {A.apply.r :f :a :b {A.new}}}}

'{def A.norms
 {def A.norms.r
  {lambda {:a :b}
   {if {A.empty? :a}
    then :b
    else {A.norms.r {A.rest :a} {A.addlast! {C.mod {A.first :a}} :b}}}}}
 {lambda {:a}
  {A.norms.r :a {A.new}}}}    
-> {def A.norms
 {def A.norms.r
  {lambda {:a :b}
   {if {A.empty? :a}
    then :b
    else {A.norms.r {A.rest :a} {A.addlast! {C.mod {A.first :a}} :b}}}}}
 {lambda {:a}
  {A.norms.r :a {A.new}}}}

'{def A.harmonics
 {def A.harmonics.r
  {lambda {:l :k :max :n}
   {if {> :k :n}
    then 
    else {if {= {A.first :l} 0}
          then
          else [f=:k,a=1/{round {/ :max {A.first :l}}}]}
         {A.harmonics.r {A.rest :l} {+ :k 1} :max :n}}}}
 {lambda {:l :N}
         {A.harmonics.r :l 0 {A.first {A.sort! > {A.duplicate :l}}} :N}}}
-> {def A.harmonics
 {def A.harmonics.r
  {lambda {:l :k :max :n}
   {if {> :k :n}
    then 
    else {if {= {A.first :l} 0} then else [f=:k,a=1/{round {/ :max {A.first :l}}}]}
         {A.harmonics.r {A.rest :l} {+ :k 1} :max :n}}}}
 {lambda {:l :N}
         {A.harmonics.r :l 0 {A.first {A.sort! > {A.duplicate :l}}} :N}}}

}

_p Note that these algorithms use complex numbers and the [[lib_complex]] library is required.

_h2 2) application

_p Find harmonics of the square curve sampled with 32 levels +1 and 32 levels -1:

_h4 1) built the sample:

{prewrap
'{def sample {A.new 
  {S.map {lambda {_} {C.new  1 0}} {S.serie 1 32}}
  {S.map {lambda {_} {C.new -1 0}} {S.serie 1 32}} }}
-> {def sample {A.new 
  {S.map {lambda {_} {C.new  1 0}} {S.serie 1 32}}
  {S.map {lambda {_} {C.new -1 0}} {S.serie 1 32}} }}   
}

_h4 2) compute the FFT:
{prewrap
'{def X {A.fft 1 {sample}}}
-> {def X {A.fft 1 {sample}}}    
}

_h4 3) show first harmonics:
{prewrap
'{A.harmonics {A.norms {X}} 11} 
-> 
{A.harmonics {A.norms {X}} 11}
}

_p We have found the five first odd sines of the curve defined by:  
{div {@ style="font:italic 1.5em courier; text-align:center;"} C(t) = Σ{sub k=0}{sup ∞} 1/k*sin(π/180*k*t) }

_p Let's plot this curve:
{pre
'{def CURVE 
 {lambda {:k :t}
  {- :t 180}
  {* 130 {+ {S.map {{lambda {:t :k}
                            {* {/ 1 :k} {sin {* {/ {PI} 180} :k :t}}} } :t} 
                   {S.serie 1 :k 2}}} }}}
-> {def CURVE 
 {lambda {:k :t} 
  {- :t 180}
  {* 127 {+ {S.map {{lambda {:t :k} {* {/ 1 :k} {sin {* {/ {PI} 180} :k :t}}} } :t} 
                   {S.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 100 0 100 0 -100 180 -100 180 0"
                                                          {stroke #888 5"}}}
  {polyline {@ points="{S.map {CURVE 1} {S.serie 0 360 2}}" {stroke #222 1}}}
  {polyline {@ points="{S.map {CURVE 3} {S.serie 0 360 2}}" {stroke #444 1}}}
  {polyline {@ points="{S.map {CURVE 5} {S.serie 0 360 2}}" {stroke #444 1}}}
  {polyline {@ points="{S.map {CURVE 7} {S.serie 0 360 2}}" {stroke #444 1}}}
  {polyline {@ points="{S.map {CURVE 9} {S.serie 0 360 2}}" {stroke #f00 3}}} 
}}} ->
{svg {@ width="600px" height="300px" style="background:#eee"}
 {g {AXES 600 300}
  {polyline {@ points="-180 0 -180 100 0 100 0 -100 180 -100 180 0" {stroke #888 5"}}}
  {polyline {@ points="{S.map {CURVE 1} {S.serie 0 360 2}}" {stroke #222 1}}}
  {polyline {@ points="{S.map {CURVE 3} {S.serie 0 360 2}}" {stroke #444 1}}}
  {polyline {@ points="{S.map {CURVE 5} {S.serie 0 360 2}}" {stroke #444 1}}}
  {polyline {@ points="{S.map {CURVE 7} {S.serie 0 360 2}}" {stroke #444 1}}}
  {polyline {@ points="{S.map {CURVE 9} {S.serie 0 360 2}}" {stroke #f00 3}}} 
}}

_ul [[http://www.rosettacode.org/wiki/Fast_Fourier_transform#lambdatalk|http://www.rosettacode.org/wiki/Fast_Fourier_transform#lambdatalk]]

_p {i Alain Marty 2020/02/02, last edit 2020/03/26}

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

pre { box-shadow:0 0 8px #000; padding:5px; }
}