Discrete Fourier transform
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:
- Implement the discrete fourier transform
- Implement the inverse fourier transform
- (optional) implement a cleaning mechanism to remove small errors introduced by floating point representation.
- 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 of length is given by:
The inverse transform is given by:
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]
Phix
include complex.e function dft(sequence x) integer N = length(x) sequence y = repeat(0,N) for k=1 to N do complex yk = complex_new(0,0) for n=1 to N do complex t = complex_new(0,-2*PI*(k-1)*(n-1)/N) yk = complex_add(yk,complex_mul(x[n],complex_exp(t))) end for y[k] = yk end for return y end function function idft(sequence y) integer N = length(y) sequence x = repeat(0,N) for n=1 to N do object xn = complex_new(0,0) for k=1 to N do complex t = complex_new(0,2*PI*(k-1)*(n-1)/N) xn = complex_add(xn,complex_mul(y[k],complex_exp(t))) end for xn = complex_div(xn,N) // clean xn to remove very small imaginary values, and round reals to 14dp if abs(complex_imag(xn))<1e-14 then xn = round(complex_real(xn),1e14) end if x[n] = xn end for return x end function sequence x = {2, 3, 5, 7, 11}, y = dft(x), z = idft(y) printf(1,"Original sequence: %v\n",{x}) printf(1,"Discrete Fourier Transform: %v\n",{apply(y,complex_sprint)}) printf(1,"Inverse Discrete Fourier Transform: %v\n",{z})
- Output:
Original sequence: {2,3,5,7,11} Discrete Fourier Transform: {"28","-3.38197+8.78402i","-5.61803+2.80017i","-5.61803-2.80017i","-3.38197-8.78402i"} Inverse Discrete Fourier Transform: {2,3,5,7,11}
Wren
<lang ecmascript>import "/complex" for Complex
var dft = Fn.new { |x|
var N = x.count var y = List.filled(N, null) for (k in 0...N) { y[k] = Complex.zero for (n in 0...N) { var t = Complex.imagMinusOne * Complex.two * Complex.pi * k * n / N y[k] = y[k] + x[n] * t.exp } } return y
}
var idft = Fn.new { |y|
var N = y.count var x = List.filled(N, null) for (n in 0...N) { x[n] = Complex.zero for (k in 0...N) { var t = Complex.imagOne * Complex.two * Complex.pi * k * n / N x[n] = x[n] + y[k] * t.exp } x[n] = x[n] / N // clean x[n] to remove very small imaginary values if (x[n].imag.abs < 1e-14) x[n] = Complex.new(x[n].real, 0) } return x
}
var x = [2, 3, 5, 7, 11] System.print("Original sequence: %(x)") for (i in 0...x.count) x[i] = Complex.new(x[i]) var y = dft.call(x) Complex.showAsReal = true // don't display the imaginary part if it's 0 System.print("\nAfter applying the Discrete Fourier Transform:") System.print(y) System.print("\nAfter applying the Inverse Discrete Fourier Transform to the above transform:") System.print(idft.call(y))</lang>
- Output:
Original sequence: [2, 3, 5, 7, 11] After applying the Discrete Fourier Transform: [28, -3.3819660112501 + 8.7840226349462i, -5.6180339887499 + 2.8001689857495i, -5.6180339887499 - 2.8001689857495i, -3.3819660112501 - 8.7840226349462i] After applying the Inverse Discrete Fourier Transform to the above transform: [2, 3, 5, 7, 11]