Julia Community 🟣

Vinod V
Vinod V

Posted on

Understanding Cooley-Tukey FFT Matrix Factorization with Julia

#Ref : https://www.cs.cornell.edu/~bindel/class/cs5220-s10/slides/FFT.pdf
#page 17
#run(`$(Base.julia_cmd())`)

using LinearAlgebra

function getBitReversed(N)
    m = Int64(log2(N))
    return [parse(Int64,reverse(bitstring(i)[end-m+1:end]),base = 2) for i in 0:N-1]
end

function getMatrices(N)
#N = 4

    t = Int64(log2(N))
    As = []
    for q = 1:t
        L = 2^q
        Lby2 = 2^(q-1)
        r = Int64(N/L)
        I_r = Matrix(1.0I,r,r)
        tw = exp.(-2im*pi*(0:Lby2-1)/L)
        D = diagm(tw)
        B_r = [Matrix(1.0I,Lby2,Lby2) D;Matrix(1.0I,Lby2,Lby2) -D]
        A_q = kron(I_r,B_r)
        push!(As,A_q)
    end
    return As
end

function test()

    N = 16

    P = 1im*zeros(N,N)

    ind = getBitReversed(N) .+ 1

    for i in Base.OneTo(N)
        P[i,ind[i]] = 1.0
    end

    As = getMatrices(N)

    F_ct = copy(P)

    for i in Base.OneTo(length(As))
        #println("A",i)
        #display(As[i])
        F_ct = As[i]*F_ct
    end

    F_N = exp.(-2im*pi*(0:N-1)*(0:N-1)'/N)

    #=
    println("F_n")
    display(F_N)

    println("F_ct")
    display(F_ct)
    =#

    println(isapprox(F_N,F_ct))
    return F_N,As,P
end

F_N,As,P = test();
Enter fullscreen mode Exit fullscreen mode

Top comments (3)

Collapse
 
oscardssmith profile image
Oscar Smith

exp.(-2im*pi*(0:Lby2-1)/L) could be cispi.(-2*(0:Lby2-1)/L) which will be slightly faster and more accurate.

Collapse
 
ohanian profile image
Steven Siew

It would have been better if you wrote some text explaining what this is all about. We CANNOT read your mind and we don't know what context or usefulness of the code.

Collapse
 
vinodv profile image
Vinod V

I have mentioned the page #13 of the reference. Anyway I will try to add some comments.