Discrete Fourier transform

Discrete Fourier transform is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The discrete Fourier transform is a linear, invertible transformation which transforms an arbitrary sequence of complex numbers to another sequence of complex numbers of the same length. The Fast Fourier transform (FFT) is an efficient implementation of this mechanism, but one which only works for sequences which have a length which is a power of 2.

The discrete Fourier transform is a useful testing mechanism to verify the correctness of code bases which use or implement the FFT.

For this task:

1. Implement the discrete fourier transform
2. Implement the inverse fourier transform
3. (optional) implement a cleaning mechanism to remove small errors introduced by floating point representation.
4. Verify the correctness of your implementation using a small sequence of integers, such as 2 3 5 7 11

The fourier transform of a sequence $x_{n}$ of length $N$ is given by: $X_{n}=\sum _{k=0}^{N-1}x_{k}\cdot e^{-i{\frac {2\pi }{N}}kn}$

The inverse transform is given by: $x_{n}={\frac {1}{N}}\sum _{k=0}^{N-1}X_{k}\cdot e^{i{\frac {2\pi }{N}}kn}$

J

Implementation: <lang j>fourier=: ] +/@:* ^@(0j_2p1 * */~@i.@# % #) ifourier=: # %~ ] +/@:* ^@(0j2p1 * */~@i.@# % #)

require'general/misc/numeric' clean=: 1e_9&round&.+.</lang>

Example use:

<lang j> clean ifourier fourier 2 3 5 7 11 2 3 5 7 11

  clean ifourier 2 * fourier 2 3 5 7 11

4 6 10 14 22

  clean ifourier 2 + fourier 2 3 5 7 11

4 3 5 7 11</lang>

Julia

<lang julia>function dft(A::AbstractArray{T,N}) where {T,N}

   F = zeros(complex(float(T)), size(A)...)
   for k in CartesianIndices(F), n in CartesianIndices(A)
       F[k] += cispi(-2 * sum(d -> (k[d] - 1) * (n[d] - 1) /
           real(eltype(F))(size(A, d)), ntuple(identity, Val{N}()))) * A[n]
   end
   return F

end

function idft(A::AbstractArray{T,N}) where {T,N}

   F = zeros(complex(float(T)), size(A)...)
   for k in CartesianIndices(F), n in CartesianIndices(A)
       F[k] += cispi(2 * sum(d -> (k[d] - 1) * (n[d] - 1) /
           real(eltype(F))(size(A, d)), ntuple(identity, Val{N}()))) * A[n]
   end
   return F ./ length(A)

end

const seq = [2, 3, 5, 7, 11]

const fseq = dft(seq)

const newseq = idft(fseq)

println("$seq =>\n$fseq =>\n$newseq =>\n", Int.(round.(newseq)))

</lang>

Output:

[2, 3, 5, 7, 11] =>
ComplexF64[28.0 + 0.0im, -3.3819660112501033 + 8.784022634946172im, -5.618033988749888 + 2.800168985749483im, -5.618033988749888 - 2.800168985749483im, -3.381966011250112 - 8.78402263494618im] =>
ComplexF64[2.0000000000000013 - 1.4210854715202005e-15im, 2.999999999999996 + 7.993605777301127e-16im, 5.000000000000002 + 2.1316282072803005e-15im, 6.999999999999998 - 8.881784197001252e-16im, 11.0 + 0.0im] =>
[2, 3, 5, 7, 11]