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

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.

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.

Vinod V

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