Discrete Fourier transform: Difference between revisions
m
→{{header|Wren}}: Minor tidy
m (comments) |
m (→{{header|Wren}}: Minor tidy) |
||
(15 intermediate revisions by 9 users not shown) | |||
Line 11:
# 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 <math>
<math>X_{n} = \sum_{k=0}^{N-1} x_{k}\cdot e^{-i \frac{2 \pi}{N} k n}</math>
The inverse transform is given by:
<math>
=={{header|Go}}==
{{trans|Wren}}
<syntaxhighlight lang="go">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()
}</syntaxhighlight>
{{out}}
<pre>
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]
</pre>
=={{header|J}}==
Implementation:
<
ifourier=: # %~ ] +/@:* ^@(0j2p1 * */~@i.@# % #)
require'general/misc/numeric'
clean=: 1e_9&round&.+.</
Example use:
<
5.6 _0.676393j1.7568 _1.12361j0.560034 _1.12361j_0.560034 _0.676393j_1.7568
clean ifourier fourier 2 3 5 7 11
Line 36 ⟶ 107:
4 6 10 14 22
clean ifourier 2 + fourier 2 3 5 7 11
4 3 5 7 11</
=={{header|jq}}==
'''Adapted from [[#Wren|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]].
<syntaxhighlight lang="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</syntaxhighlight>
{{out}}
<pre>
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]
</pre>
=={{header|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
<
F = zeros(complex(float(T)), size(A)...)
for k in CartesianIndices(F), n in CartesianIndices(A)
Line 58 ⟶ 189:
end
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)
println("$seq2 =>\n$fseq =>\n$newseq =>\n", Int.(round.(newseq)))
</syntaxhighlight>{{out}}
<pre>
[2, 3, 5, 7, 11] =>
Line 71 ⟶ 204:
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]
</pre>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">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</syntaxhighlight>
{{out}}
<pre>{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.}</pre>
=={{header|Nim}}==
{{trans|Wren}}
<syntaxhighlight lang="nim">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(", ")</syntaxhighlight>
{{out}}
<pre>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</pre>
=={{header|Perl}}==
<syntaxhighlight lang="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;
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|Phix}}==
{{trans|Wren}}
<!--<
<span style="color: #008080;">include</span> <span style="color: #004080;">complex</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>
Line 115 ⟶ 357:
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Discrete Fourier Transform: %v\n"</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">apply</span><span style="color: #0000FF;">(</span><span style="color: #000000;">y</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">complex_sprint</span><span style="color: #0000FF;">)})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Inverse Discrete Fourier Transform: %v\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">z</span><span style="color: #0000FF;">})</span>
<!--</
{{out}}
<pre>
Line 121 ⟶ 363:
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}
</pre>
=={{header|Python}}==
<syntaxhighlight lang="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 ))
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
Line 134 ⟶ 476:
* <tt>$n.round($r)</tt> returns ''$n'' rounded to the nearest ''$r''. <tt>(13/16).round(1/4)</tt> is <tt>3/4</tt> .
* .narrow changes a Numeric's type to a "narrower" type, when no precision would be lost.<br><tt>8/2</tt> is <tt>4</tt>, but is type Rat (rational). <tt>(8/2).narrow</tt> is also <tt>4</tt>, but type Int.
<syntaxhighlight lang="raku"
my @exp := ( $pos_neg_i * tau / +@x * $k ) «*« @x.keys;
return sum @x »*« 𝑒 «**« @exp;
Line 152 ⟶ 494:
say '';
warn "Round-trip failed" unless ( clean(@x) Z== clean(@x_idft) ).all;
}</
{{out}}
<pre>
Line 164 ⟶ 506:
=={{header|Wren}}==
{{libheader|Wren-complex}}
<
var dft = Fn.new { |x|
Line 203 ⟶ 545:
System.print(y)
System.print("\nAfter applying the Inverse Discrete Fourier Transform to the above transform:")
System.print(idft.call(y))</
{{out}}
|