2016 Update
This commit is contained in:
Tina Müller
parent
948b86eafa
commit
dcf5d15da3
7965 changed files with 139854 additions and 31002 deletions
|
|
@ -1,10 +1,13 @@
|
|||
{{omit from|GUISS}}
|
||||
The purpose of this task is to calculate the FFT (Fast Fourier Transform)
|
||||
of an input sequence. <br>
|
||||
|
||||
;Task:
|
||||
Calculate the FFT (<u>F</u>ast <u>F</u>ourier <u>T</u>ransform) of an input sequence.
|
||||
|
||||
The most general case allows for complex numbers at the input
|
||||
and results in a sequence of equal length, again of complex numbers.
|
||||
If you need to restrict yourself to real numbers, the output should
|
||||
be the magnitude (i.e. sqrt(re²+im²)) of the complex result.
|
||||
|
||||
The classic version is the recursive Cooley–Tukey FFT. [http://en.wikipedia.org/wiki/Cooley–Tukey_FFT_algorithm Wikipedia] has pseudocode for that.
|
||||
The classic version is the recursive Cooley–Tukey FFT. [http://en.wikipedia.org/wiki/Cooley–Tukey_FFT_algorithm Wikipedia] has pseudo-code for that.
|
||||
Further optimizations are possible but not required.
|
||||
<br><br>
|
||||
|
|
|
|||
|
|
@ -27,23 +27,24 @@ void fft(cplx buf[], int n)
|
|||
_fft(buf, out, n, 1);
|
||||
}
|
||||
|
||||
|
||||
void show(const char * s, cplx buf[]) {
|
||||
printf("%s", s);
|
||||
for (int i = 0; i < 8; i++)
|
||||
if (!cimag(buf[i]))
|
||||
printf("%g ", creal(buf[i]));
|
||||
else
|
||||
printf("(%g, %g) ", creal(buf[i]), cimag(buf[i]));
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
PI = atan2(1, 1) * 4;
|
||||
cplx buf[] = {1, 1, 1, 1, 0, 0, 0, 0};
|
||||
|
||||
void show(const char * s) {
|
||||
printf("%s", s);
|
||||
for (int i = 0; i < 8; i++)
|
||||
if (!cimag(buf[i]))
|
||||
printf("%g ", creal(buf[i]));
|
||||
else
|
||||
printf("(%g, %g) ", creal(buf[i]), cimag(buf[i]));
|
||||
}
|
||||
|
||||
show("Data: ");
|
||||
show("Data: ", buf);
|
||||
fft(buf, 8);
|
||||
show("\nFFT : ");
|
||||
show("\nFFT : ", buf);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,2 +1,41 @@
|
|||
Data: 1 1 1 1 0 0 0 0
|
||||
FFT : 4 (1, -2.41421) 0 (1, -0.414214) 0 (1, 0.414214) 0 (1, 2.41421)
|
||||
#include <stdio.h>
|
||||
#include <Accelerate/Accelerate.h>
|
||||
|
||||
void fft(DSPComplex buf[], int n) {
|
||||
float inputMemory[2*n];
|
||||
float outputMemory[2*n];
|
||||
// half for real and half for complex
|
||||
DSPSplitComplex inputSplit = {inputMemory, inputMemory + n};
|
||||
DSPSplitComplex outputSplit = {outputMemory, outputMemory + n};
|
||||
|
||||
vDSP_ctoz(buf, 2, &inputSplit, 1, n);
|
||||
|
||||
vDSP_DFT_Setup setup = vDSP_DFT_zop_CreateSetup(NULL, n, vDSP_DFT_FORWARD);
|
||||
|
||||
vDSP_DFT_Execute(setup,
|
||||
inputSplit.realp, inputSplit.imagp,
|
||||
outputSplit.realp, outputSplit.imagp);
|
||||
|
||||
vDSP_ztoc(&outputSplit, 1, buf, 2, n);
|
||||
}
|
||||
|
||||
|
||||
void show(const char *s, DSPComplex buf[], int n) {
|
||||
printf("%s", s);
|
||||
for (int i = 0; i < n; i++)
|
||||
if (!buf[i].imag)
|
||||
printf("%g ", buf[i].real);
|
||||
else
|
||||
printf("(%g, %g) ", buf[i].real, buf[i].imag);
|
||||
printf("\n");
|
||||
}
|
||||
|
||||
int main() {
|
||||
DSPComplex buf[] = {{1,0}, {1,0}, {1,0}, {1,0}, {0,0}, {0,0}, {0,0}, {0,0}};
|
||||
|
||||
show("Data: ", buf, 8);
|
||||
fft(buf, 8);
|
||||
show("FFT : ", buf, 8);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,50 +0,0 @@
|
|||
#include <stdio.h>
|
||||
#include <math.h>
|
||||
#include <complex.h>
|
||||
|
||||
double PI;
|
||||
typedef double complex cplx;
|
||||
|
||||
void _fft(cplx buf[], cplx out[], int n, int step)
|
||||
{
|
||||
if (step < n) {
|
||||
_fft(out, buf, n, step * 2);
|
||||
_fft(out + step, buf + step, n, step * 2);
|
||||
|
||||
for (int i = 0; i < n; i += 2 * step) {
|
||||
cplx t = cexp(-I * PI * i / n) * out[i + step];
|
||||
buf[i / 2] = out[i] + t;
|
||||
buf[(i + n)/2] = out[i] - t;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void fft(cplx buf[], int n)
|
||||
{
|
||||
cplx out[n];
|
||||
for (int i = 0; i < n; i++) out[i] = buf[i];
|
||||
|
||||
_fft(buf, out, n, 1);
|
||||
}
|
||||
|
||||
|
||||
void show(const char * s, cplx buf[]) {
|
||||
printf("%s", s);
|
||||
for (int i = 0; i < 8; i++)
|
||||
if (!cimag(buf[i]))
|
||||
printf("%g ", creal(buf[i]));
|
||||
else
|
||||
printf("(%g, %g) ", creal(buf[i]), cimag(buf[i]));
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
PI = atan2(1, 1) * 4;
|
||||
cplx buf[] = {1, 1, 1, 1, 0, 0, 0, 0};
|
||||
|
||||
show("Data: ", buf);
|
||||
fft(buf, 8);
|
||||
show("\nFFT : ", buf);
|
||||
|
||||
return 0;
|
||||
}
|
||||
95
Task/Fast-Fourier-transform/Java/fast-fourier-transform.java
Normal file
95
Task/Fast-Fourier-transform/Java/fast-fourier-transform.java
Normal file
|
|
@ -0,0 +1,95 @@
|
|||
import static java.lang.Math.*;
|
||||
|
||||
public class FastFourierTransform {
|
||||
|
||||
public static int bitReverse(int n, int bits) {
|
||||
int reversedN = n;
|
||||
int count = bits - 1;
|
||||
|
||||
n >>= 1;
|
||||
while (n > 0) {
|
||||
reversedN = (reversedN << 1) | (n & 1);
|
||||
count--;
|
||||
n >>= 1;
|
||||
}
|
||||
|
||||
return ((reversedN << count) & ((1 << bits) - 1));
|
||||
}
|
||||
|
||||
static void fft(Complex[] buffer) {
|
||||
|
||||
int bits = (int) (log(buffer.length) / log(2));
|
||||
for (int j = 1; j < buffer.length / 2; j++) {
|
||||
|
||||
int swapPos = bitReverse(j, bits);
|
||||
Complex temp = buffer[j];
|
||||
buffer[j] = buffer[swapPos];
|
||||
buffer[swapPos] = temp;
|
||||
}
|
||||
|
||||
for (int N = 2; N <= buffer.length; N <<= 1) {
|
||||
for (int i = 0; i < buffer.length; i += N) {
|
||||
for (int k = 0; k < N / 2; k++) {
|
||||
|
||||
int evenIndex = i + k;
|
||||
int oddIndex = i + k + (N / 2);
|
||||
Complex even = buffer[evenIndex];
|
||||
Complex odd = buffer[oddIndex];
|
||||
|
||||
double term = (-2 * PI * k) / (double) N;
|
||||
Complex exp = (new Complex(cos(term), sin(term)).mult(odd));
|
||||
|
||||
buffer[evenIndex] = even.add(exp);
|
||||
buffer[oddIndex] = even.sub(exp);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
public static void main(String[] args) {
|
||||
double[] input = {1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0};
|
||||
|
||||
Complex[] cinput = new Complex[input.length];
|
||||
for (int i = 0; i < input.length; i++)
|
||||
cinput[i] = new Complex(input[i], 0.0);
|
||||
|
||||
fft(cinput);
|
||||
|
||||
System.out.println("Results:");
|
||||
for (Complex c : cinput) {
|
||||
System.out.println(c);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
class Complex {
|
||||
public final double re;
|
||||
public final double im;
|
||||
|
||||
public Complex() {
|
||||
this(0, 0);
|
||||
}
|
||||
|
||||
public Complex(double r, double i) {
|
||||
re = r;
|
||||
im = i;
|
||||
}
|
||||
|
||||
public Complex add(Complex b) {
|
||||
return new Complex(this.re + b.re, this.im + b.im);
|
||||
}
|
||||
|
||||
public Complex sub(Complex b) {
|
||||
return new Complex(this.re - b.re, this.im - b.im);
|
||||
}
|
||||
|
||||
public Complex mult(Complex b) {
|
||||
return new Complex(this.re * b.re - this.im * b.im,
|
||||
this.re * b.im + this.im * b.re);
|
||||
}
|
||||
|
||||
@Override
|
||||
public String toString() {
|
||||
return String.format("(%f,%f)", re, im);
|
||||
}
|
||||
}
|
||||
67
Task/Fast-Fourier-transform/Lua/fast-fourier-transform.lua
Normal file
67
Task/Fast-Fourier-transform/Lua/fast-fourier-transform.lua
Normal file
|
|
@ -0,0 +1,67 @@
|
|||
-- operations on complex number
|
||||
complex = {__mt={} }
|
||||
|
||||
function complex.new (r, i)
|
||||
local new={r=r, i=i or 0}
|
||||
setmetatable(new,complex.__mt)
|
||||
return new
|
||||
end
|
||||
|
||||
function complex.__mt.__add (c1, c2)
|
||||
return complex.new(c1.r + c2.r, c1.i + c2.i)
|
||||
end
|
||||
|
||||
function complex.__mt.__sub (c1, c2)
|
||||
return complex.new(c1.r - c2.r, c1.i - c2.i)
|
||||
end
|
||||
|
||||
function complex.__mt.__mul (c1, c2)
|
||||
return complex.new(c1.r*c2.r - c1.i*c2.i,
|
||||
c1.r*c2.i + c1.i*c2.r)
|
||||
end
|
||||
|
||||
function complex.expi (i)
|
||||
return complex.new(math.cos(i),math.sin(i))
|
||||
end
|
||||
|
||||
function complex.__mt.__tostring(c)
|
||||
return "("..c.r..","..c.i..")"
|
||||
end
|
||||
|
||||
|
||||
-- Cooley–Tukey FFT (in-place, divide-and-conquer)
|
||||
-- Higher memory requirements and redundancy although more intuitive
|
||||
function fft(vect)
|
||||
local n=#vect
|
||||
if n<=1 then return vect end
|
||||
-- divide
|
||||
local odd,even={},{}
|
||||
for i=1,n,2 do
|
||||
odd[#odd+1]=vect[i]
|
||||
even[#even+1]=vect[i+1]
|
||||
end
|
||||
-- conquer
|
||||
fft(even);
|
||||
fft(odd);
|
||||
-- combine
|
||||
for k=1,n/2 do
|
||||
local t=even[k] * complex.expi(-2*math.pi*(k-1)/n)
|
||||
vect[k] = odd[k] + t;
|
||||
vect[k+n/2] = odd[k] - t;
|
||||
end
|
||||
return vect
|
||||
end
|
||||
|
||||
function toComplex(vectr)
|
||||
vect={}
|
||||
for i,r in ipairs(vectr) do
|
||||
vect[i]=complex.new(r)
|
||||
end
|
||||
return vect
|
||||
end
|
||||
|
||||
-- test
|
||||
data = toComplex{1, 1, 1, 1, 0, 0, 0, 0};
|
||||
|
||||
print("orig:", unpack(data))
|
||||
print("fft:", unpack(fft(data)))
|
||||
|
|
@ -2,13 +2,13 @@ sub fft {
|
|||
return @_ if @_ == 1;
|
||||
my @evn = fft( @_[0, 2 ... *] );
|
||||
my @odd = fft( @_[1, 3 ... *] ) Z*
|
||||
map &cis, (0, 2 * pi / @_ ... *);
|
||||
map &cis, (0, tau / @_ ... *);
|
||||
return flat @evn »+« @odd, @evn »-« @odd;
|
||||
}
|
||||
|
||||
my @seq = ^16;
|
||||
my $cycles = 3;
|
||||
my @wave = map { sin( 2*pi * $_ / @seq * $cycles ) }, @seq;
|
||||
my @wave = map { sin( tau * $_ / @seq * $cycles ) }, @seq;
|
||||
say "wave: ", @wave.fmt("%7.3f");
|
||||
|
||||
say "fft: ", fft(@wave)».abs.fmt("%7.3f");
|
||||
|
|
|
|||
|
|
@ -1,70 +1,69 @@
|
|||
/*REXX pgm performs a fast Fourier transform (FFT) on a set of complex numbers*/
|
||||
numeric digits length( pi() ) - 1 /*limited by the PI function result. */
|
||||
arg data /*ARG verb uppercases the DATA from CL.*/
|
||||
if data='' then data=1 1 1 1 0 /*Not specified? Then use the default.*/
|
||||
size=words(data); pad=left('',6) /*PAD: for indenting and padding SAYs.*/
|
||||
do p=0 until 2**p>=size ; end /*number of args exactly a power of 2? */
|
||||
do j=size+1 to 2**p;data=data 0; end /*add zeroes to DATA 'til a power of 2.*/
|
||||
size=words(data); ph=p%2; call hdr /*╔═════════════════════════════╗*/
|
||||
/* [↓] TRANSLATE allows I&J*/ /*║ Numbers in data can be in ║*/
|
||||
do j=0 for size /*║ seven formats: real ║*/
|
||||
_=translate(word(data,j+1), 'J', "I") /*║ real,imag ║*/
|
||||
parse var _ #.1.j '' $ 1 ',' #.2.j /*║ ,imag ║*/
|
||||
if $=='J' then parse var #.1.j #2.j , /*║ nnnJ ║*/
|
||||
"J" #.1.j /*║ nnnj ║*/
|
||||
do m=1 for 2; #.m.j=word(#.m.j 0,1) /*║ nnnI ║*/
|
||||
end /*m*/ /* [↑] ommited part?*/ /*║ nnni ║*/
|
||||
/*╚═════════════════════════════╝*/
|
||||
say pad " FFT in " center(j+1,7) pad fmt(#.1.j) fmt(#.2.j,'i')
|
||||
end /*j*/
|
||||
/*REXX program performs a fast Fourier transform (FFT) on a set of complex numbers. */
|
||||
numeric digits length( pi() ) - 1 /*limited by the PI function result. */
|
||||
arg data /*ARG verb uppercases the DATA from CL.*/
|
||||
if data='' then data=1 1 1 1 0 /*Not specified? Then use the default.*/
|
||||
size=words(data); pad=left('',6) /*PAD: for indenting and padding SAYs.*/
|
||||
do p=0 until 2**p>=size ; end /*number of args exactly a power of 2? */
|
||||
do j=size+1 to 2**p;data=data 0; end /*add zeroes to DATA 'til a power of 2.*/
|
||||
size=words(data); ph=p%2; call hdr /*╔═══════════════════════════╗*/
|
||||
/* [↓] TRANSLATE allows I & J*/ /*║ Numbers in data can be in ║*/
|
||||
do j=0 for size /*║ seven formats: real ║*/
|
||||
_=translate( word(data,j+1), 'J', "I") /*║ real,imag ║*/
|
||||
parse var _ #.1.j '' $ 1 "," #.2.j /*║ ,imag ║*/
|
||||
if $=='J' then parse var #.1.j #2.j "J" #.1.j /*║ nnnJ ║*/
|
||||
/*║ nnnj ║*/
|
||||
do m=1 for 2; #.m.j= word(#.m.j 0, 1) /*║ nnnI ║*/
|
||||
end /*m*/ /*omitted part? [↑] */ /*║ nnni ║*/
|
||||
/*╚═══════════════════════════╝*/
|
||||
say pad ' FFT in ' center(j+1, 7) pad fmt(#.1.j) fmt(#.2.j, "i")
|
||||
end /*j*/
|
||||
say
|
||||
tran=pi()*2/2**p; !.=0; hp=2**p%2; A=2**(p-ph); ptr=A; dbl=1
|
||||
tran=pi()*2 / 2**p; !.=0; hp=2**p %2; A=2**(p-ph); ptr=A; dbl=1
|
||||
say
|
||||
do p-ph; halfPtr=ptr%2
|
||||
do i=halfPtr by ptr to A-halfPtr; _=i-halfPtr; !.i=!._+dbl
|
||||
end /*i*/
|
||||
dbl=dbl*2; ptr=halfPtr
|
||||
end /*p-ph*/
|
||||
do p-ph; halfPtr=ptr%2
|
||||
do i=halfPtr by ptr to A-halfPtr; _=i-halfPtr; !.i=!._+dbl
|
||||
end /*i*/
|
||||
dbl=dbl*2; ptr=halfPtr
|
||||
end /*p-ph*/
|
||||
|
||||
do j=0 to 2**p%4; cmp.j=cos(j*tran); _=hp - j; cmp._= -cmp.j
|
||||
_=hp + j; cmp._= -cmp.j
|
||||
end /*j*/
|
||||
do j=0 to 2**p%4; cmp.j=cos(j*tran); _=hp - j; cmp._= -cmp.j
|
||||
_=hp + j; cmp._= -cmp.j
|
||||
end /*j*/
|
||||
B=2**ph
|
||||
|
||||
do i=0 for A; q=i * B
|
||||
do j=0 for B; h=q+j; _=!.j*B+!.i; if _<=h then iterate
|
||||
parse value #.1._ #.1.h #.2._ #.2.h with #.1.h #.1._ #.2.h #.2._
|
||||
end /*j*/ /* [↑] swap two sets of values.*/
|
||||
end /*i*/
|
||||
|
||||
dbl=1; do p ; w=hp % dbl
|
||||
do k=0 for dbl ; Lb=w * k ; Lh=Lb + 2**p % 4
|
||||
do j=0 for w ; a=j * dbl * 2 + k ; b= a + dbl
|
||||
r=#.1.a; i=#.2.a ; c1=cmp.Lb * #.1.b ; c4=cmp.Lb * #.2.b
|
||||
c2=cmp.Lh * #.2.b ; c3=cmp.Lh * #.1.b
|
||||
#.1.a=r + c1 - c2 ; #.2.a=i + c3 + c4
|
||||
#.1.b=r - c1 + c2 ; #.2.b=i - c3 - c4
|
||||
end /*j*/
|
||||
end /*k*/
|
||||
dbl=dbl+dbl
|
||||
end /*p*/
|
||||
do i=0 for A; q=i * B
|
||||
do j=0 for B; h=q+j; _=!.j*B+!.i; if _<=h then iterate
|
||||
parse value #.1._ #.1.h #.2._ #.2.h with #.1.h #.1._ #.2.h #.2._
|
||||
end /*j*/ /* [↑] swap two sets of values. */
|
||||
end /*i*/
|
||||
dbl=1
|
||||
do p ; w=hp % dbl
|
||||
do k=0 for dbl ; Lb=w * k ; Lh=Lb + 2**p % 4
|
||||
do j=0 for w ; a=j * dbl * 2 + k ; b= a + dbl
|
||||
r=#.1.a; i=#.2.a ; c1=cmp.Lb * #.1.b ; c4=cmp.Lb * #.2.b
|
||||
c2=cmp.Lh * #.2.b ; c3=cmp.Lh * #.1.b
|
||||
#.1.a=r + c1 - c2 ; #.2.a=i + c3 + c4
|
||||
#.1.b=r - c1 + c2 ; #.2.b=i - c3 - c4
|
||||
end /*j*/
|
||||
end /*k*/
|
||||
dbl=dbl+dbl
|
||||
end /*p*/
|
||||
call hdr
|
||||
do i=0 for size
|
||||
say pad " FFT out " center(i+1,7) pad fmt(#.1.i) fmt(#.2.i,'j')
|
||||
end /*i*/ /*numbers are shown with 10 digs [↑] */
|
||||
exit /*stick a fork in it, we're all done. */
|
||||
/*────────────────────────────────────────────────────────────────────────────*/
|
||||
cos: procedure; parse arg x; q=r2r(x)**2; z=1; _=1; p=1
|
||||
do k=2 by 2; _=-_*q/(k*(k-1)); z=z+_; if z=p then leave; p=z; end; return z
|
||||
/*────────────────────────────────────────────────────────────────────────────*/
|
||||
fmt: procedure; parse arg y,j; y=y/1 /*transforms complex numbers for looks.*/
|
||||
if abs(y)<'1e-'digits()%4 then y=0; if y=0 & j\=='' then return ''
|
||||
y=format(y,,10); if pos(.,y)\==0 then y=strip(y,'T',0)
|
||||
y=strip(y,,.); if y>=0 then y=' 'y; return left(y||j, 12)
|
||||
/*────────────────────────────────────────────────────────────────────────────*/
|
||||
hdr: _='───data─── num real─part imaginary─part'; say pad _
|
||||
say pad translate(_, " "copies('═',256), " "xrange()); return
|
||||
/*────────────────────────────────────────────────────────────────────────────*/
|
||||
pi: return 3.141592653589793238462643383279502884197169399375105820974944592308
|
||||
/*────────────────────────────────────────────────────────────────────────────*/
|
||||
r2r: return arg(1) // (pi()*2) /*reduce the radians to a unit circle. */
|
||||
end /*i*/ /*[↑] #s are shown with 10 decimal digs*/
|
||||
exit /*stick a fork in it, we're all done. */
|
||||
/*──────────────────────────────────────────────────────────────────────────────────────*/
|
||||
cos: procedure; parse arg x; q=r2r(x)**2; z=1; _=1; p=1 /*bare bones COS. */
|
||||
do k=2 by 2; _=-_*q/(k*(k-1)); z=z+_; if z=p then leave; p=z; end; return z
|
||||
/*──────────────────────────────────────────────────────────────────────────────────────*/
|
||||
fmt: procedure; parse arg y,j; y=y/1 /*prettifies complex numbers for show. */
|
||||
if abs(y) < '1e-'digits()%4 then y=0; if y=0 & j\=='' then return ''
|
||||
y=format(y, , 10); if pos(.,y)\==0 then y=strip(y, 'T', 0)
|
||||
y=strip(y, , .); if y>=0 then y=' 'y; return left(y || j, 12)
|
||||
/*──────────────────────────────────────────────────────────────────────────────────────*/
|
||||
hdr: _='───data─── num real─part imaginary─part'; say pad _
|
||||
say pad translate(_, " "copies('═',256), " "xrange()); return
|
||||
/*──────────────────────────────────────────────────────────────────────────────────────*/
|
||||
pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164062862
|
||||
/*──────────────────────────────────────────────────────────────────────────────────────*/
|
||||
r2r: return arg(1) // ( pi()*2 ) /*reduce the radians to a unit circle. */
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue