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();