Let’s implement the FFT algorithm in julia. We need a bit reversal function, and we create a single-pass function that can be iterated. These are two utility functions for bit reversing integers and arrays indexed by integers.

""" note pad=log2(n)"""
bitrev=function(num,n,pad)
foo=[2^(k-1) for k in 1:pad]
tmp=reverse(digits(num,base=2,pad=pad))
sum(foo .* tmp)
end

bitrevV=function(f)
m=length(f)
pad=Int(log2(m))
indices=[bitrev(j,m,pad) for j in 0:(m-1)]
[f[m+1] for m in indices]
end

Next we create the actual FFT function that iterates a pass function

"""N must be a power of two"""
sfft=function(f,L)
N=length(f)
powr=Int64(log2(N))-1
frev=bitrevV(f)
K=[j-N/2 for j in 0:N-1]
g=[0.0+0.0*im for j in 0:N-1]
for j in 0:(N-1)
ω=exp(-2π*K[j+1]*im/N)
z=Complex.(frev)
for a in powr:-1:0
    z=pass2(z,ω,a)
end
g[j+1]=z[1]
end
K,g
end


pass2=function(f,ω,exponent)
n=length(f)
factor=ω^(2^exponent)
retval=[f[2*m+1]+factor*f[2*m+2] for m in 0:(Int64(n/2)-1)]
return(retval)
end

Try it out on some data, we sample a superposition of cosines

using Gaston
"""make a signal"""
NN=2^(8)
t = [2π*j/(NN-1) for j in 0:NN-1]
signal = cos.(3.0 .* t) .+ 0.5*cos.(9.0 .* t)

x,y=sfft(signal,5)
plot(x,abs.(y),w=:lp)

We expect delta functions, we get delta functions…

Next time: using the julia FFTW library.