Fourier transform?

Hi guys.
I uninstalled CODEA some months ago because I hadn’t upgraded to IOS5 and it was just too frustating to read about the awsome new features and not being able to test them… Now I have updated and it is just too cool to see all the way you’ve done in just a few months. Cargo bot is amazing! Thanks a lot for all that good work.
Now, since you can do so many miracles, i wondered if you would some time implement 2d Fourier transform in an complex numbers extention? It is not very difficult, and then i coud use codea to implement a series of cool image processing tools, as a showcase of recent image processing progresses. Image processing is my business and i’d like to port it on the ipad. I could use Javascript but LUA+CODEA is so much more efficient for making nice GUIs, and also you provide simple APIs to output result images which are not available in Javascript. I think your platform is the only one on Ipad with such completeness (hope Big Brother is looking elsewhere…).
Anyway, thanks again for all the pleasure you give to us casual programmers, and keep your enthousiasm.
Best regards

Thanks for the kind words @Jmv38

Are you thinking of something like Lapack?

Perhaps we could build an API on top of Apple’s Accelerate framework and our image type. I haven’t touched this sort of image processing for a while, it might be an area @Dylan can contribute towards.

@Simeon: thanks for answering. Lapack seems to be a linear algebra extension oriented on matrix calculations. I thought of something much simpler with: complex numbers, -,/,+,*,= operations, modulus() and arg() fonctions, and 1d complex FFT on rows or column of a matrix (image). Some morphological functions might be necesary too. Maybe i should program the functions i need in lua, regardless on the performance, an then ask if some of you guys would be interested in making a CODEA library out of it to accelerate the output?

I’ve written a library for complex numbers with the standard operations, not FFT as yet but wouldn’t be hard to add. You can get it from http://www.math.ntnu.no/~stacey/documents/Codea/ComplexNumbers/ (though that might be an old version - I’ll check later today).

@andrew_stacey looks great! I could start from that. Thanks!

Just checked - that’s the same as my current version.

I have some vector and matrix code as well, but that’s not online yet.

I am currently working on a FFT algorithm. I can’t promiss it 'll be fast, but will share the code.
It could use some optimising, so fire away. I’m always eager to learn from my mistakes



function FastFourierTransform(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 inversetransform(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 = FastFourierTransform(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

PS How do I pass the values of an array to a function so that the original isn’t affected?
In the implementation above the function works directly on the array passed to it.
That’s memory efficient, but perhaps not always wanted

Maybe try make a shallow copy:

dataCopy = {}
for k, v in ipairs(data) do
    dataCopy[k] = i
end

@PTN thanks for sharing. I haven’t carefully studied your code, but it seems to be a direct calculation of each coef - the hard way in term of operation amount (O(n2)). I would inplement the butterfly algorithm which is much faster (O(n.log(n))). I’ve also read the MIT is about to release an even faster algo, but no news about them yet.

This version does use the butterfly scheme to reduce the number of calculations. it could be made faster by stopping earlier (when the data has four points).
Another way is by assuming the imaginary part of the input to be zero and then make use of the fact that the frequency spectrum is symmetrical.

It seems to me that these inner calculations resemble vector operations, so perhaps using vectors may be faster.
I will try to implement these points and post them

With copying your data, you’ll have to do a deep copy because everything in lua is a pointer. Thus:

~% lua
Lua 5.1.4  Copyright (C) 1994-2008 Lua.org, PUC-Rio
> a = {3,4}
> b = a
> b[1] = 5
> print(a)
table: 0x7fca1bc0b9c0
> print(unpack(a))
5	4
> 

So since your input is a table of tables, a shallow copy will still point to the same tables. The easiest way would be when doing the bit reversal: have a new table hold the sorted values rather than reusing the old one.

Also, if you’re going to use symbolic names for the real and imaginary parts, please don’t use r and i! With complex numbers, r is used for the length and i is just a square root of -1. Use x and y!

@PTN sorry, i hadn’t noticed the butterfly. Please don’t assume the imaginary part is 0! If the input image is real, the processed FFT itsefl is not, so the inverse FFT operation has to work with complex numbers => there is only little gain to speed up the direct FFT only… Thanks anyway.

I happened to have some code for manipulating vectors, complex numbers, and binary numbers lying around so I put them together and added an FFT routine. It’s the “standard” one and extends the input vector to one of length 2^n padding with zeros.

Available as part of my general library project.

@Andrew_Stacey thanks a lot! Where exactly can i find your ‘general library project’?

I’ve just posted about it. The short version is: http://www.math.ntnu.no/~stacey/code/CodeaLibrary

@Andrew_Stacey Thanks. I see the fft is in the ‘‘vector.lua’’ file.

Great stuff with the library and all @Andrew_Stacey. Soooo much to learn from that…
Optimizing for real-valued inputs is giving me a headache. The speed of the FFT as-is is fast enough for me for now. Perhaps when applying it in 2D speed wil be an issue again.
The r/i to x/y naming conversion is a good point. Thanks !