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 $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]

Phix

Translation of: Wren

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}

Raku

This task could be done with a loop of maps like
@X[k] = sum @x.kv.map: -> \n, \v { v * exp( -i * tau / N * k * n ) }
, but it is a better fit for Raku's concurrent Hyper-operators.
In the DFT formula, N becomes +@x, the element count. We would usually omit the plus sign and get the same result from numeric context, but that gets confusing to read when mixed with hyper-ops like »*« .
The exponents of DFT and IDFT only differ from each other in the sign of i, and the calculation of any given elements only differ by the k index of the element being output, so the inner loop is cleanly shareable between DFT and IDFT.
Euler's constant, imaginary unit, pi, and 2pi are built-in, with ASCII and Unicode spellings: e 𝑒 i pi π tau τ .
$one_thing «op« @many_things vectorizes the op operation: 5 «+« (1,2,3) is (6,7,8) .
@many_things »op« @many_more_things distributes the op operation pairwise: (5,6,7) »+« (1,2,3) is (6,8,10) .
(list)».method applies the method to each element of the list.
@array.keys, when @array has e.g. 4 elements, is 0,1,2,3 .
$n.round($r) returns $n rounded to the nearest $r. (13/16).round(1/4) is 3/4 .
.narrow changes a Numeric's type to a "narrower" type, when no precision would be lost.
8/2 is 4, but is type Rat (rational). (8/2).narrow is also 4, but type Int.

<lang perl6>sub ft_inner ( @x, $k, $pos_neg_i where * == i|-i ) {

   my @exp := ( $pos_neg_i * tau / +@x * $k ) «*« @x.keys;
   return sum @x »*« 𝑒 «**« @exp;

} sub dft ( @x ) { return @x.keys.map: { ft_inner( @x, $_, -i ) } } sub idft ( @x ) { return @x.keys.map: { ft_inner( @x, $_, i ) / +@x } } sub clean ( @x ) { @x».round(1e-12)».narrow }

my @tests = ( 1, 2-i, -i, -1+2i ),

           ( 2,   3,  5,     7, 11 ),

for @tests -> @x {

   my @x_dft  =  dft(@x);
   my @x_idft = idft(@x_dft);

   say .key.fmt('%6s:'), .value.&clean.fmt('%5s', ', ') for :@x, :@x_dft, :@x_idft;
   say ;
   warn "Round-trip failed" unless ( clean(@x) Z== clean(@x_idft) ).all;

}</lang>

Output:

     x:    1,  2-1i,  0-1i, -1+2i
 x_dft:    2, -2-2i,  0-2i,  4+4i
x_idft:    1,  2-1i,  0-1i, -1+2i

     x:    2,     3,     5,     7,    11
 x_dft:   28, -3.38196601125+8.784022634946i, -5.61803398875+2.800168985749i, -5.61803398875-2.800168985749i, -3.38196601125-8.784022634946i
x_idft:    2,     3,     5,     7,    11

Wren

Library: Wren-complex

<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]