Skip to content

Instantly share code, notes, and snippets.

Created Oct 7, 2014
Embed
What would you like to do?
Gists Codea Upload
--# fft
-- fft
function realToCplx(iTable)
cplx = {}
for i,v in ipairs(iTable) do
cplx[i-1] = { r=v, i=0}
end
return cplx
end
function getField(iTable,field)
t = {}
for i=0,#iTable do
t[i+1] = iTable[i][field]
end
return t
end
function ri2mp(iTable)
t = {}
for k=0,#iTable do
local c = iTable[k]
local r,i = c.r,c.i
local out = {}
out.m = math.sqrt( r*r + i*i)
if out.m < 1/10000000000 then
out.p = 0
else
out.p = math.atan2(i,r)
end
t[k] = out
end
return t
end
function fft(data)
-- data contains n points of complex numbers
-- data[i].r is the real part
-- data[i].i is the imaginary part
-- n must be a power of 2
-- the array you pass will contain the transformed data after this function call
-- the original data will be lost. Be sure to store the data elsewhere if you need it again.
-- If the data size is the wrong length it will return false
local N = #data + 1
if N == 0 or N == nil then
print("no data")
return false
end
local m = 0
local e = 0
m,e = math.frexp(N)
if m ~= 0.5 then
print("size of data not a power of 2")
return false
end
local PI = math.pi
local NM1 = N-1
local NM2 = N-2
local ND2 = N/2
local M = e - 1
local J = ND2
local K = 0
local Tr = 0
local Ti = 0
-- do the bit reversal sorting part
for i = 1,NM2 do
if i < J then
Tr = data[J].r
Ti = data[J].i
data[J].r = data[i].r
data[J].i = data[i].i
data[i].r = Tr
data[i].i = Ti
end
K = ND2
while K <= J do
J = J - K
K = K/2
end
J = J + K
end
local LE = 0
local LE2 = 0
local UR = 0
local UI = 0
local SR = 0
local SI = 0
local JM = 0
local IP = 0
for L = 1,M do
LE = 2^L
LE2 = LE/2
UR = 1
UI = 0
SR = math.cos(PI/LE2)
SI = - math.sin(PI/LE2)
for J = 1,LE2 do
JM1 = J-1
for I = JM1, NM1, LE do
IP = I+LE2
Tr = data[IP].r * UR - data[IP].i * UI
Ti = data[IP].r * UI + data[IP].i * UR
data[IP].r = data[I].r - Tr
data[IP].i = data[I].i - Ti
data[I].r = data[I].r + Tr
data[I].i = data[I].i + Ti
end
Tr = UR
UR = Tr*SR - UI*SI
UI = Tr*SI + UI*SR
end
end
return true
end
function ifft(data)
-- data contains n points of complex numbers
-- data[i].r is the real part
-- data[i].i is the imaginary part
-- n must be a power of 2
-- the array you pass will contain the transformed data after this function call
-- the original data will be lost. Be sure to store the data elsewhere if you need it again.
-- If the data size is the wrong length it will return false.
local NN = #data + 1
for k = 0,NN-1 do
data[k].r = data[k].r
data[k].i = -data[k].i
end
local ssucces = fft(data)
for i = 0, NN-1 do
data[i].r = data[i].r / NN
data[i].i = -data[i].i / NN
end
return ssucces
end
--# Main
-- forum fft
-- Use this function to perform your initial setup
function setup()
parameter.watch( "math.floor(1/DeltaTime)" )
n = 64
parameter.integer("start",0,n-1,0)
parameter.integer("stop",0,n-1,n-1)
parameter.integer("NbCyclesPerImage",0,n/2,4)
Plot:setup(1,3)
p1 = Plot(1,1)
p1.title = "Original data"
c2 = p1:newCurve()
p2 = Plot(1,2)
p2.title = "Fourier Transform: modulus"
c3 = p2:newCurve()
p3 = Plot(1,3)
p3.title = "Fourier Transform: phase"
c4 = p3:newCurve()
update()
p1:autoScaleFromCurve(c2)
p2:autoScaleFromCurve(c2)
p2:scaleY(0,n/2)
p3:autoScaleFromCurve(c2)
p3:scaleY(-4,4)
end
function update()
local count = n / NbCyclesPerImage
local test = sinCurve1(n, count, math.pi/2 + ElapsedTime, start,stop)
local cplx = realToCplx( test )
c2:update( test )
fft(cplx)
local mp = ri2mp( cplx )
c3:update( getField(mp,"m") )
c4:update( getField(mp,"p") )
end
-- This function gets called once every frame
function draw()
background(40, 40, 50)
p1:draw()
p2:draw()
p3:draw()
update()
end
function sinCurve(n,pointsPerCycle, phi)
local y = {}
local pi = math.pi
local phi = phi or 0
for i = 0,n-1 do
y[i+1] = math.sin( pi * 2 * i / pointsPerCycle + phi)
end
return y
end
function sinCurve1(n,pointsPerCycle, phi,start,stop)
local y = {}
local pi = math.pi
local phi = phi or 0
local gate,v
for i = 0,n-1 do
v = math.sin( pi * 2 * i / pointsPerCycle + phi)
if (i<start) or (i>stop) then gate = 0 else gate = 1 end
y[i+1] = v * gate
end
return y
end
local random = math.random
function randomCurve(n)
local y = {}
local pi = math.pi
local phi = phi or 0
for i = 0,n-1 do
y[i+1] = random()*2-1
end
return y
end
--# Style
Style = class()
function Style:init()
-- stores control functions for each field that needs it
self.control = {}
end
function Style:set(field,value)
-- first check if the value is acceptable, and modify it if needed
local f = self.control[field]
if f then value = f(self,value) end
self[field] = value
end
--# Curve
Curve = class()
function Curve:init(parentPlot)
self.parentPlot = parentPlot
self.data = {}
self:defaultStyle()
end
function Curve:defaultStyle(data)
self.style = Style()
self.style.pointRadius = 15
self.style.pointFill = color(0,0)
self.style.pointStroke = color(0)
self.style.pointStrokeWidth = 1.5
self.style.lineStroke = color(255, 0, 0, 255)
self.style.lineStrokeWidth = 5
end
function Curve:update(data)
self.data = data
end
function Curve:draw()
local r = self.style.pointRadius
local pf = self.style.pointFill
local ps = self.style.pointStroke
local psw = self.style.pointStrokeWidth
local Ls = self.style.lineStroke
local Lsw = self.style.lineStrokeWidth
ellipseMode(CENTER)
fill(pf)
stroke(Ls)
strokeWidth(Lsw)
local sx,sy = self.parentPlot.sx, self.parentPlot.sy
local x0,y0, x1,y1
local f = self.process
for i,v in ipairs( self.data ) do
local x,y = i,v
if f then x,y = f(x,y) end
x1,y1 = x*sx , y*sy
if i>1 then
line(x0,y0,x1,y1)
end
x0,y0 = x1,y1
end
stroke(ps)
strokeWidth(psw)
for i,v in ipairs( self.data ) do
local x,y = i,v
if f then x,y = f(x,y) end
x1,y1 = x*sx , y*sy
ellipse( x1,y1 , r,r)
end
end
function Curve:drawTitle(x,y,w)
if not self.title then return end
fill(self.style.lineStroke )
textWrapWidth(w)
text(self.title, x,y)
end
function Curve:touched(touch)
-- Codea does not automatically call this method
end
--# Plot
Plot = class()
-- number of plots on the screen
function Plot:setup(Nx,Ny)
Plot.Nx = Nx or 1
Plot.Ny = Ny or 1
end
function Plot:init(nx,ny)
self:defaultStyle()
if Plot.Nx == nil then Plot:setup() end
local margin = self.style.margin
local wplot, hplot = (WIDTH-margin)/Plot.Nx, (HEIGHT-margin)/Plot.Ny
self.x = wplot * (Plot.Nx- nx) + margin
self.w = wplot - margin
self.y = hplot * (Plot.Ny - ny) + margin
self.h = hplot - margin
self.curves = {}
self.ox = self.w/2
self.oy = self.h/2
self.sx = 10
self.sy = 10
end
function Plot:defaultStyle()
self.style = Style()
self.style.margin = 2
self.style.fill = color(220, 255)
self.style.stroke = color(0)
self.style.strokeWidth = 2
self.style.titleFont = "Arial-ItalicMT"
self.style.titleFontSize = 25
self.style.titleFontColor = color(0)
self.style.titleOffsetY = 0
self.style.titleOffsetX = 0
self.style.legend = false
self.style.legendWidth = 200
end
Plot.defaultLineStrokes = {
color(255, 0, 0, 255),
color(0, 0, 255, 255),
color(21, 178, 21, 255),
color(0, 255),
color(255, 163, 0, 255),
color(0, 217, 255, 255)
}
function Plot:newCurve()
local curve = Curve(self)
table.insert( self.curves, curve)
local n = #self.curves
local c = Plot.defaultLineStrokes[n]
if c then curve.style:set("lineStroke",c) end
return curve
end
function Plot:autoScaleFromCurve(curve)
local data = curve.data
if #data == 0 then return end
local mini,maxi = 1,#data
self:scaleX(mini,maxi)
local mini,maxi = data[1],data[1]
for i,y in ipairs(data) do
if y < mini then mini = y end
if y > maxi then maxi = y end
end
local d = (maxi - mini) /10
self:scaleY(mini-d,maxi+d)
end
function Plot:scaleX(mini,maxi)
local d = (maxi-mini)
if d > 0 then
local w = 0
if self.style.legend then
w = self.style.legendWidth
self.style.titleOffsetX = -w/2
end
self.sx = (self.w - w -10 )/ d
end
self.ox = - mini * self.sx
end
function Plot:scaleY(mini,maxi)
local d = (maxi-mini)
if d > 0 then
local h = 0
if self.title then h = (self.style.titleFontSize + self.style.titleOffsetY) *1.5 end
self.sy = (self.h - h ) / d
end
self.oy = - mini * self.sy
end
function Plot:draw()
resetMatrix()
resetStyle()
-- plot decorations
translate(self.x,self.y)
local s = self.style
-- plot frame
fill( s.fill )
stroke( s.stroke )
strokeWidth( s.strokeWidth )
rect(0,0,self.w,self.h)
-- plot title
if self.title then
fill( s.titleFontColor )
font( s.titleFont )
local h = s.titleFontSize
fontSize(h)
textMode(CENTER)
text( self.title, self.w/2+ s.titleOffsetX, self.h - h + s.titleOffsetY)
end
-- curves
for i,curve in pairs(self.curves) do
pushMatrix()
translate(self.ox,self.oy)
curve:draw()
popMatrix()
if self.style.legend then
local w = self.style.legendWidth
local x = self.w - w/2
local dh = self.h/(#self.curves+1)
local y = dh * i
curve:drawTitle( x,y,w)
end
end
end
function Plot:touched(t)
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment