Gists Codea Upload

Project.lua

--# 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