Fast Fourier transform: Difference between revisions
(J) |
|||
Line 1: | Line 1: | ||
{{draft task}} |
{{draft task}} |
||
The purpose of this task is to calculate the FFT (Fast Fourier Transform) 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. Further optimizations are possible but not required. |
The purpose of this task is to calculate the FFT (Fast Fourier Transform) 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. Further optimizations are possible but not required. |
||
=={{header|J}}== |
|||
From [[j:Essays/FFT]]: |
|||
<lang j>cube =: ($~ q:@#) :. , |
|||
roots =: +@]^:(0>[) (_1^2%]) ^ i.@-: NB. The elegant version, but less accurate |
|||
rou =: [: (, j.) [: (, (j.~%:0.5) , |."1&.:+.@|.@}.) [: ^@o.@:j. i.@(%&8) % -: NB. roots of unity |
|||
floop =: 4 : 'for_r. i.#$x do. (y=.{."1 y) ] x=.(+/x) ,&,:"r (-/x)*y end.' |
|||
fft =: ( ] floop&.cube rou@#) f. :. ifft |
|||
ifft =: (# %~ ] floop&.cube +@rou@#) f. :. fft</lang> |
|||
Example: |
|||
<lang j>fmt =: [:, sprintf~&'%7.3f '@| |
|||
fmt"1 (,:fft) 1 o. 3r8p1 * i.16 |
|||
0.000 0.924 0.707 0.383 1.000 0.383 0.707 0.924 0.000 0.924 0.707 0.383 1.000 0.383 0.707 0.924 |
|||
0.000 0.000 0.000 8.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 8.000 0.000 0.000</lang> |
|||
Note that fmt only displays the magnitude of the result (and shows only 3 decimals of precision). |
|||
=={{header|Perl 6}}== |
=={{header|Perl 6}}== |
Revision as of 16:32, 8 February 2011
The purpose of this task is to calculate the FFT (Fast Fourier Transform) 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. Wikipedia has pseudocode for that. Further optimizations are possible but not required.
J
From j:Essays/FFT:
<lang j>cube =: ($~ q:@#) :. , roots =: +@]^:(0>[) (_1^2%]) ^ i.@-: NB. The elegant version, but less accurate rou =: [: (, j.) [: (, (j.~%:0.5) , |."1&.:+.@|.@}.) [: ^@o.@:j. i.@(%&8) % -: NB. roots of unity floop =: 4 : 'for_r. i.#$x do. (y=.{."1 y) ] x=.(+/x) ,&,:"r (-/x)*y end.' fft =: ( ] floop&.cube rou@#) f. :. ifft ifft =: (# %~ ] floop&.cube +@rou@#) f. :. fft</lang>
Example:
<lang j>fmt =: [:, sprintf~&'%7.3f '@|
fmt"1 (,:fft) 1 o. 3r8p1 * i.16 0.000 0.924 0.707 0.383 1.000 0.383 0.707 0.924 0.000 0.924 0.707 0.383 1.000 0.383 0.707 0.924 0.000 0.000 0.000 8.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 8.000 0.000 0.000</lang>
Note that fmt only displays the magnitude of the result (and shows only 3 decimals of precision).
Perl 6
<lang perl6>sub fft {
return @_ if @_ == 1; my @evn = fft( @_[0,2...^* >= @_] ); my @odd = fft( @_[1,3...^* >= @_] ); my $twd = 2i * pi / @_; # twiddle factor @odd »*=« (^@odd).map(* * $twd)».exp; return @evn »+« @odd, @evn »-« @odd;
}
my @seq = ^16; my $cycles = 3; my @wave = (@seq »*» (2*pi / @seq * $cycles))».sin; say "wave: ", @wave.fmt("%7.3f");
say "fft: ", fft(@wave)».abs.fmt("%7.3f");</lang>
Output:
wave: 0.000 0.924 0.707 -0.383 -1.000 -0.383 0.707 0.924 0.000 -0.924 -0.707 0.383 1.000 0.383 -0.707 -0.924 fft: 0.000 0.000 0.000 8.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 8.000 0.000 0.000