Testing the code for the DFT algorithm written in julia
In dft1, I naively translate the C algorithm into julia without take advantage of julia's complex values
The length in dft1 is twice the length that it should be becuase the complex values are being stored as indata[n] is real and indata[n+1] is imaginary, where \(n = 0,2,4,...\)
I implemented my algrorithm as dft2 which was optimized to use the complex values of julia ( also because this made it easier to implement ). I used linspace function to generate a 1000 points of the Sine function and ran it through the algorithm. dft2 generated questionable out put. Comparing this to Jim's version I implemented in julia, dft3, it is clear the outout of my dft2 is not correct. I'll need to go through it again and look at what I did wrong. It may be an error in translating it from C to julia but it could also be an error in my math.
function dft1(indata, outdata)
sumreal = sumimag = 0
N = length(indata)/2.0
for i in [0:N-1]
for k in [0:N-1]
sumreal += indata[2*k] * cos(2 * pi * k * i / N) + indata[2*k+1] * sin(2 * pi * k * i / N)
sumimag += -indata[2*k] * sin(2 * pi * k * i / N) + indata[2*k+1] * cos(2 * pi * k * i / N)
end
outdata[2*i] = sumreal
outdata[2*i+1] = sumimag
end
end
function dft2(indata)
sumreal = sumimag = 0
N = length(indata)
outdata = Array(Complex, N)
for i in [1:N]
for k in [1:N]
sumreal += real(indata[k]) * cos(2 * pi * k * i / N) + imag(indata[k]) * sin(2 * pi * k * i / N)
sumimag += -real(indata[k]) * sin(2 * pi * k * i / N) + imag(indata[k]) * cos(2 * pi * k * i / N)
end
outdata[i] = complex(sumreal,sumimag)
end
return outdata
end
Z = dft2(Y)
function dft3(indata)
N = length(indata)
omega = exp( -1im * 2pi / N)
outdata = Array(Complex, N)
for i in [1:N]
for k in [1:N]
outdata[i] = indata[i] * omega^(i*k)
end
end
return outdata
end
B = dft3(Y)
using PyPlot
plot(real(B),imag(B))
title("Test plot of Sine")
xlabel("Frequency")
Well, that didn't work :
plot([1:1000]/1000.0,real(B))
title("Test plot of Sine in Frequency domain 2")
xlabel("Frequency")
wait, what? Now, I'm kinda confused. Why did I just get the sin wave back? I thought there would be a spike at one frequency...
plot([1:1000]/1000.0,imag(B))
title("Test plot of Sine in frequency domain 3")
xlabel("Frequency")
Okay, so the imaginary plot did not reval anything. It looks like garbage due to arithmetic artifacts. In the Test plot of Sine in frequency domain 2 I think what its telling me is that the sine wave is being approximated by a series of different sine waves which happend to have thier frequencies defined by a sine function. Unless, I did something wrong and that is total bullshit.
Z
for i in [1:length(y)]
y[i] = complex(y[i],0)
end
Y = Array(Complex, length(y))
for i in [1:length(y)]
Y[i] = complex(y[i],0)
end
Y
x = linspace(0, 2pi, 1000)
y = sin(x)
f = g = [1:10]
f
g
x = 5
[1:x-1]
real(complex(1,-1))
complex(1,-1)
real(-1im)
real(-1im)
reset
real
real(-1im)
real(-1im)
foo = Complex[1:10]
foo = Array(Complex, 10)
exp(1)
exp( -1im )
exp( -1im * 2pi )