Discrete Fourier transform

From Rosetta Code
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 of length is given by:


The inverse transform is given by:

Go

Translation of: Wren
package main

import (
    "fmt"
    "math"
    "math/cmplx"
)

func dft(x []complex128) []complex128 {
    N := len(x)
    y := make([]complex128, N)
    for k := 0; k < N; k++ {
        for n := 0; n < N; n++ {
            t := -1i * 2 * complex(math.Pi*float64(k)*float64(n)/float64(N), 0)
            y[k] += x[n] * cmplx.Exp(t)
        }
    }
    return y
}

func idft(y []complex128) []float64 {
    N := len(y)
    x := make([]complex128, N)
    for n := 0; n < N; n++ {
        for k := 0; k < N; k++ {
            t := 1i * 2 * complex(math.Pi*float64(k)*float64(n)/float64(N), 0)
            x[n] += y[k] * cmplx.Exp(t)
        }
        x[n] /= complex(float64(N), 0)
        // clean x[n] to remove very small imaginary values
        if math.Abs(imag(x[n])) < 1e-14 {
            x[n] = complex(real(x[n]), 0)
        }
    }
    z := make([]float64, N)
    for i, c := range x {
        z[i] = float64(real(c))
    }
    return z
}

func main() {
    z := []float64{2, 3, 5, 7, 11}
    x := make([]complex128, len(z))
    fmt.Println("Original sequence:", z)
    for i, n := range z {
        x[i] = complex(n, 0)
    }
    y := dft(x)
    fmt.Println("\nAfter applying the Discrete Fourier Transform:")
    fmt.Printf("%0.14g", y)
    fmt.Println("\n\nAfter applying the Inverse Discrete Fourier Transform to the above transform:")
    z = idft(y)
    fmt.Printf("%0.14g", z)
    fmt.Println()
}
Output:
Original sequence: [2 3 5 7 11]

After applying the Discrete Fourier Transform:
[(28+0i) (-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]

J

Implementation:

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

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

Example use:

   clean fourier 2 3 5 7 11
5.6 _0.676393j1.7568 _1.12361j0.560034 _1.12361j_0.560034 _0.676393j_1.7568
   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

jq

Adapted from Wren

Works with: jq

Works with gojq, the Go implementation of jq

This entry uses a "complex" module, such as is available at Arithmetic/Complex#jq.

include "complex"; # a reminder
 
def dft:
  . as $x
  | length as $N
  | reduce range (0; $N) as $k ([]; # y
      .[$k] = [0,0] # Complex.zero
      | reduce range( 0; $N) as $n (.;
          ([[0, -1], [2,0], [pi,0], $k, $n, invert($N) ] | multiply) as $t
          | .[$k] = plus(.[$k]; multiply($x[$n]; exp($t))) ) ) ;

# Input: an array of Complex
def idft:
  . as $y
  | length as $N
  | reduce range(0; $N) as $n ([];
      .[$n] = [0,0]
      | reduce range(0; $N) as $k (.;
          ( [ 2, pi, [0,1], $k, $n, invert($N)] | multiply) as $t
          | .[$n] = plus(.[$n];  multiply($y[$k]; exp($t))) )
        | .[$n] = divide(.[$n]; $N) );


def task:
  def abs: if . < 0 then -. else . end;
  # round, and remove very small imaginary values altogether
  def tidy:
    (.[0] | round) as $round
    | if (.[0]| (. - $round) | abs < 1e-14) then .[0] = $round else . end
    | if .[1]|abs < 1e-14 then .[0] else . end;

  [2, 3, 5, 7, 11] as $x
  | "Original sequence: \($x)",
    (reduce range(0; $x|length) as $i ( []; .[$i] = [$x[$i], 0])
     | dft
     | "\nAfter applying the Discrete Fourier Transform:",
       .,
       "\nAfter applying the Inverse Discrete Fourier Transform to this value: \(
          idft | map(tidy))" )
     ;

task
Output:
Original sequence: [2,3,5,7,11]

After applying the Discrete Fourier Transform:
[[28,0],[-3.3819660112501078,8.784022634946174],[-5.618033988749895,2.8001689857494747],[-5.618033988749892,-2.8001689857494823],[-3.381966011250096,-8.784022634946178]]

After applying the Inverse Discrete Fourier Transform to this value: [2,3,5,7,11]


Julia

The cispi function was added in Julia 1.6. The cispi of x is e to the power of (pi times i times x). Other syntax, such as {T,N} in the function signature and real() and ntuple() in the loop, is designed to support arbitrary dimensional arrays and to cue the compiler so as to have mostly constants in the loop at runtime, speeding run time with large arrays.

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

seq = [2, 3, 5, 7, 11]
fseq = dft(seq)
newseq = idft(fseq)
println("$seq =>\n$fseq =>\n$newseq =>\n", Int.(round.(newseq)))

seq2 = [2 7; 3 11; 5 13]
fseq = dft(seq2)
newseq = idft(fseq)
println("$seq2 =>\n$fseq =>\n$newseq =>\n", Int.(round.(newseq)))
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]
[2 7; 3 11; 5 13] =>
ComplexF64[41.0 + 0.0im -21.0 + 0.0im; -7.000000000000002 + 3.46410161513775im 3.0000000000000053 + 7.105427357601002e-15im; -6.9999999999999964 - 3.4641016151377606im 3.0000000000000053 + 7.105427357601002e-15im] =>
ComplexF64[2.000000000000002 + 5.921189464667501e-16im 6.999999999999997 - 4.144832625267251e-15im; 2.9999999999999996 - 8.881784197001252e-16im 11.000000000000002 + 1.0362081563168128e-15im; 4.999999999999999 + 0.0im 13.000000000000002 + 1.924386576016938e-15im] =>
[2 7; 3 11; 5 13]

Mathematica/Wolfram Language

ClearAll[DFT, IDFT]
DFT[x_List] := Module[{length},
  length = Length[x];
  N@Table[Sum[x[[n + 1]] Exp[-I 2 Pi k n/length], {n, 0, length - 1}], {k, 0, length - 1}]
 ]
IDFT[X_List] := Module[{length},
  length = Length[X];
  N@Table[Sum[X[[k + 1]] Exp[-I 2 Pi k n/length], {k, 0, length - 1}]/length, {n, 0, length - 1}]
 ]
DFT[{2, 3, 5, 7, 11}]
IDFT[%] // Chop
Output:
{28., -3.38197 + 8.78402 I, -5.61803 + 2.80017 I, -5.61803 - 2.80017 I, -3.38197 - 8.78402 I}
{2., 11., 7., 5., 3.}

Nim

Translation of: Wren
import complex, math, sequtils, strutils

func dft(x: openArray[Complex64]): seq[Complex64] =
  let N = x.len
  result.setLen(N)
  for k in 0..<N:
    for n in 0..<N:
      let t = complex64(0, -2 * PI * float(k) * float(n) / float(N))
      result[k] += x[n] * exp(t)


func idft(y: openArray[Complex64]): seq[Complex64] =
  let N = y.len
  result.setLen(N)
  let d = complex64(float(N))
  for n in 0..<N:
    for k in 0..<N:
      let t = complex64(0, 2 * PI * float(k) * float(n) / float(N))
      result[n] += y[k] * exp(t)
    result[n] /= d
    # Clean result[n] to remove very small imaginary values.
    if abs(result[n].im) < 1e-14: result[n].im = 0.0


func `$`(c: Complex64): string =
  result = c.re.formatFloat(ffDecimal, precision = 2)
  if c.im != 0:
    result.add if c.im > 0: "+" else: ""
    result.add c.im.formatFloat(ffDecimal, precision = 2) & 'i'


when isMainModule:

  let x = [float 2, 3, 5, 7, 11].mapIt(complex64(it))
  echo "Original sequence: ", x.join(", ")
  let y = dft(x)
  echo "Discrete Fourier transform: ", y.join(", ")
  echo "Inverse Discrete Fourier Transform: ", idft(y).join(", ")
Output:
Original sequence: 2.00, 3.00, 5.00, 7.00, 11.00
Discrete Fourier transform: 28.00, -3.38+8.78i, -5.62+2.80i, -5.62-2.80i, -3.38-8.78i
Inverse Discrete Fourier Transform: 2.00, 3.00, 5.00, 7.00, 11.00

Perl

# 20210616 Perl programming solution

use strict;
use warnings;
use feature 'say';

use Math::Complex;
use constant  PI => 4 * atan2(1, 1);
use constant ESP =>  1e10; # net worth, the imaginary part
use constant EPS => 1e-10; # the reality part

sub dft {
   my $n = scalar ( my @in = @_ );
   return map {
      my $s=0;
      for my $k (0 .. $n-1) { $s += $in[$k] * exp(-2*i * PI * $k * $_ / $n) }
      $_ = $s;
   } (0 .. $n-1);
}

sub idft {
   my $n = scalar ( my @in = @_ );
   return map {
      my $s=0;
      for my $k (0 .. $n-1) { $s += $in[$k] * exp(2*i * PI * $k * $_ / $n) }
      my $t = $s/$n;
      $_ = abs(Im $t) < EPS ? Re($t) : $t 
   } (0 .. $n-1);
}

say 'Original sequence                  : ', join ', ', my @series = ( 2, 3, 5, 7, 11 );
say 'Discrete Fourier transform         : ', join ', ', my @dft = dft @series;
say 'Inverse Discrete Fourier Transform : ', join ', ', idft @dft;
Output:
Original sequence                  : 2, 3, 5, 7, 11
Discrete Fourier transform         : 28, -3.38196601125011+8.78402263494617i,-5.6180339887499+2.80016898574947i, -5.61803398874989-2.80016898574948i, -3.3819660112501-8.78402263494618i
Inverse Discrete Fourier Transform : 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}

Python

# 20220918 Python programming solution
import cmath


def dft( x ):
    """Takes N either real or complex signal samples, yields complex DFT bins.  Assuming the input
    waveform is filtered to only contain signals in the bandwidth B range -B/2:+B/2 around baseband
    frequency MID, and is frequency shifted (divided by) your baseband frequency MID, and is sampled
    at the Nyquist rate R: given N samples, the result contains N signal frequency component bins:

        index:      0                           N/2                   N-1
        baseband:  [MID+] [MID+] ... [MID+]    [MID+/-]  [MID+]  ... [MID+]
        frequency:  DC     1B/N   (N/2-1)B/N   (N/2)B/N  (1-N/2)B/N  -1B/N

    """
    N                           = len( x )
    result                      = []
    for k in range( N ):
        r                       = 0
        for n in range( N ):
            t                   = -2j * cmath.pi * k * n / N
            r                  += x[n] * cmath.exp( t )
        result.append( r )
    return result


def idft( y ):
    """Inverse DFT on complex frequency bins."""
    N                           = len( y )
    result                      = []
    for n in range( N ):
        r                       = 0
        for k in range( N ):
            t                   = 2j * cmath.pi * k * n / N
            r                  += y[k] * cmath.exp( t )
        r                      /= N+0j
        result.append( r )
    return result


if __name__ == "__main__":
    x                           = [ 2, 3, 5, 7, 11 ]
    print( "vals:   " + ' '.join( f"{f:11.2f}" for f in x ))
    y                           = dft( x )
    print( "DFT:    " + ' '.join( f"{f:11.2f}" for f in y ))
    z                           = idft( y )
    print( "inverse:" + ' '.join( f"{f:11.2f}" for f in z ))
    print( " - real:" + ' '.join( f"{f.real:11.2f}" for f in z ))

    N                           = 8
    print( f"Complex signals, 1-4 cycles in {N} samples; energy into successive DFT bins" )
    for rot in (0, 1, 2, 3, -4, -3, -2, -1):    # cycles; and bins in ascending index order
        if rot > N/2:
            print( "Signal change frequency exceeds sample rate and will result in artifacts")
        sig                     = [
            # unit-magnitude complex samples, rotated through 2Pi 'rot' times, in N steps
            cmath.rect(
                1, cmath.pi*2*rot/N*i
            )
            for i in range( N )
        ]
        print( f"{rot:2} cycle" + ' '.join( f"{f:11.2f}" for f in sig ))
        dft_sig                 = dft( sig )
        print( f"  DFT:  " + ' '.join( f"{f:11.2f}" for f in dft_sig ))
        print( f"   ABS: " + ' '.join( f"{abs(f):11.2f}" for f in dft_sig ))
Output:
vals:          2.00        3.00        5.00        7.00       11.00
DFT:    28.00+0.00j -3.38+8.78j -5.62+2.80j -5.62-2.80j -3.38-8.78j
inverse: 2.00-0.00j  3.00-0.00j  5.00+0.00j  7.00+0.00j 11.00+0.00j
 - real:       2.00        3.00        5.00        7.00       11.00
Complex signals, 1-4 cycles in 8 samples; energy into successive DFT bins
 0 cycle 1.00+0.00j  1.00+0.00j  1.00+0.00j  1.00+0.00j  1.00+0.00j  1.00+0.00j  1.00+0.00j  1.00+0.00j
  DFT:   8.00+0.00j -0.00+0.00j -0.00-0.00j -0.00+0.00j  0.00-0.00j -0.00-0.00j -0.00-0.00j  0.00+0.00j
   ABS:        8.00        0.00        0.00        0.00        0.00        0.00        0.00        0.00
 1 cycle 1.00+0.00j  0.71+0.71j  0.00+1.00j -0.71+0.71j -1.00+0.00j -0.71-0.71j -0.00-1.00j  0.71-0.71j
  DFT:  -0.00-0.00j  8.00+0.00j -0.00+0.00j  0.00-0.00j -0.00-0.00j  0.00+0.00j  0.00+0.00j  0.00+0.00j
   ABS:        0.00        8.00        0.00        0.00        0.00        0.00        0.00        0.00
 2 cycle 1.00+0.00j  0.00+1.00j -1.00+0.00j -0.00-1.00j  1.00-0.00j  0.00+1.00j -1.00+0.00j -0.00-1.00j
  DFT:  -0.00+0.00j -0.00-0.00j  8.00+0.00j -0.00-0.00j -0.00-0.00j  0.00-0.00j  0.00-0.00j  0.00-0.00j
   ABS:        0.00        0.00        8.00        0.00        0.00        0.00        0.00        0.00
 3 cycle 1.00+0.00j -0.71+0.71j -0.00-1.00j  0.71+0.71j -1.00+0.00j  0.71-0.71j  0.00+1.00j -0.71-0.71j
  DFT:  -0.00-0.00j  0.00+0.00j -0.00+0.00j  8.00+0.00j -0.00+0.00j -0.00-0.00j  0.00+0.00j  0.00-0.00j
   ABS:        0.00        0.00        0.00        8.00        0.00        0.00        0.00        0.00
-4 cycle 1.00-0.00j -1.00-0.00j  1.00+0.00j -1.00-0.00j  1.00+0.00j -1.00-0.00j  1.00+0.00j -1.00-0.00j
  DFT:   0.00-0.00j  0.00-0.00j  0.00-0.00j  0.00-0.00j  8.00+0.00j -0.00-0.00j  0.00-0.00j -0.00-0.00j
   ABS:        0.00        0.00        0.00        0.00        8.00        0.00        0.00        0.00
-3 cycle 1.00-0.00j -0.71-0.71j -0.00+1.00j  0.71-0.71j -1.00-0.00j  0.71+0.71j  0.00-1.00j -0.71+0.71j
  DFT:  -0.00+0.00j  0.00-0.00j  0.00+0.00j -0.00-0.00j  0.00-0.00j  8.00+0.00j -0.00-0.00j -0.00+0.00j
   ABS:        0.00        0.00        0.00        0.00        0.00        8.00        0.00        0.00
-2 cycle 1.00-0.00j  0.00-1.00j -1.00-0.00j -0.00+1.00j  1.00+0.00j  0.00-1.00j -1.00-0.00j -0.00+1.00j
  DFT:  -0.00-0.00j -0.00-0.00j  0.00-0.00j  0.00+0.00j  0.00-0.00j  0.00-0.00j  8.00+0.00j -0.00+0.00j
   ABS:        0.00        0.00        0.00        0.00        0.00        0.00        8.00        0.00
-1 cycle 1.00-0.00j  0.71-0.71j  0.00-1.00j -0.71-0.71j -1.00-0.00j -0.71+0.71j -0.00+1.00j  0.71+0.71j
  DFT:  -0.00+0.00j  0.00-0.00j -0.00-0.00j  0.00-0.00j  0.00-0.00j  0.00-0.00j  0.00-0.00j  8.00+0.00j
   ABS:        0.00        0.00        0.00        0.00        0.00        0.00        0.00        8.00

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.
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;
}
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
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))
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]