Fast Fourier transform: Difference between revisions

Content deleted Content added
Line 1,474:
<lang scheme>
 
1) the Cfft function (direct & inverse)fft
See http://lambdaway.free.fr/lambdaspeech/?view=lispology4 for more details
 
{def CUrootfft
1) the Cfft function (direct & inverse)
{lambda {:xs :f}
{if {= {list.length :f} 1}
then :f
else {combine :s {fft :s {evens :f}}
{seriefft 1:s {#.lengthodds :x} 2f}}} }}}
 
{def Cfftcombine
{lambda {:xs :sev :od}
{ifplusminus :ev {>rotate :s 0 {#list.length :xod} 1:od}}}}
then {Ccombine :x :s
{Cfft {#.evens :x} :s}
{Cfft {#.odds :x} :s} 0 {#.length :x}}
else :x}}}
 
{def Ccombineplusminus
{lambda {:xa :s :ev :od :k :Nb}
{iflist.append {<list.map Cadd :k {/a :N 2}b}
then {Ccombine {Cplusminus :x :s :ev :od :k :N} :s :ev :od {+list.map Csub :ka 1:b}}} :N}
else :x}}}
 
{def Cplusminusrotate
{lambda {:x :s :evk :odN :k :Nf}
{letif { {:x :x} {:k :k} {:Nlist.null? :Nf}
then nil
{:evk {#.get :ev :k}}
else {cons {:rootCmul {CUrootcar :s :od :k :N}f} }
{#.set! {#.set! :x {+Cexp :k{Cnew 0 {/ {* :Ns 2}{PI} {Csub :evkk} :rootN}}} }
{rotate :s {+ :k 1} :kN {Caddcdr :evk :rootf}} } }}}
 
2) four functions onfor arrayslists
{def CUroot
{lambda {:s :od :k :N}
{Cmul {#.get :od :k} {Cexp {Cnew 0 {/ {* -2 :s {PI} :k} :N}}}} }}
 
We add to the existing {lambda talk}'s list primitives a small set of functions required by the function fft.
2) four functions on arrays
 
{def #.evens
{lambda {:xl}
{if {list.null? :l}
{#.new {map {{lambda {:x :i} {#.get :x :i}} :x}
then nil
{serie 0 {- {#.length :x} 1} 2}}} }}
else {cons {car :l} {evens {cdr {cdr :l}}}}}}}
 
{def #.odds
{lambda {:xl}
{if {list.null? {cdr :l}}
{#.new {map {{lambda {:x :i} {#.get :x :i}} :x}
then nil
{serie 1 {#.length :x} 2}}} }}
else {cons {car {cdr :l}} {odds {cdr {cdr :l}}}}}}}
 
{def #list.sample2Complexmap
{def list.map.r
{lambda {:x}
{#.new {map {{lambda {:xf :i} {Cnew {#.geta :x :i} 0}}b :xc}
{serie 0if {- {#list.lengthnull? :x} 1}}}} }a}
then :c
else {list.map.r :f {cdr :a} {cdr :b}
{serie 0 {-cons {#.length:f {car :xa} 1}{car 2:b}} :c}} }}}
{lambda {:sf :oda :k :Nb}
{list.map.r :f {list.reverse :a} {list.reverse :b} nil}}}
 
{def #list.duplicateappend
{def list.append.r
{lambda {:a}
{#.slicelambda {:a 0 {#.length :a}} }b}
{if {list.null? :b}
then :a
else {list.append.r {cons {car :b} :a} {cdr :b}}}}}
{lambda {:a :b}
{list.append.r :b {list.reverse :a}} }}
 
3) a small libraryfunctions for Cnumbers
 
{lambda talk} has no primitive functions working on complex numbers. We add the minimal set required by the function fft.
 
{def Cnew
Line 1,552 ⟶ 1,564:
 
{def Cexp
{lambda {:x}
{cons {* {exp {car :x}} {cos {cdr :x}}}
{* {exp {car :x}} {sin {cdr :x}}}} }}
 
{def Clist
{lambda {:as}
{#list.new {map {{lambda {:x :i} {#.getcons :xi :i0}} :xs}}}}
 
4) testing
 
Applying the fft function on such a sample (1 1 1 1 0 0 0 0) where numbers have been promoted as complex
We define a reduced sample of the square curve [1 1 1 1 0 0 0 0], translate it into Cnumbers, then apply the FFT {Cfft x 1} and the inverse FFT {Cfft x -1} to retrieve the initial sample.
 
{def sample {#.new 1 1 1 1 0 0 0 0}}
-> [1 1 1 1 0 0 0 0]
 
{def x {#.sample2Complex {sample}}}
-> [(1 0),(1 0),(1 0),(1 0),(0 0),(0 0),(0 0),(0 0)]
 
{def Xlist.disp {Cfftfft {#.duplicate-1 {x}}Clist 1 1 1 1 0 0 0 0}}} ->
[(4 0),
(1 -2.414213562373095),
(0 0),
(1 -0.4142135623730949),
(0 0),
(0.9999999999999999 0.4142135623730949),
(0 0),
(0.9999999999999997 2.414213562373095)]
 
[(4 0),
{def x' {Cfft {#.duplicate {X}} -1}} ->
(1 -2.414213562373095),
[(8 0),
(0 0),
(8 -4.440892098500626e-16),
(1 -0.4142135623730949),
(8 1.9915985002059197e-16),
(0 0),
(8 -4.440892098500626e-16),
(0.9999999999999999 0.4142135623730949),
(4.440892098500626e-16 0),
(0 0),
(4.440892098500626e-16 4.440892098500626e-16),
(0.9999999999999997 2.414213562373095)]
(0 -1.9915985002059197e-16),
(4.440892098500626e-16 4.440892098500626e-16)]
 
SeeA more usefull example can be seen in http://lambdaway.free.fr/lambdaspeech/?view=lispology4 for more detailszorg
The last expression x' is equivalent to the initial sample where values are divided by 8 and values close to zero are displayed as 0.
 
</lang>