Fast Fourier transform

From Rosetta Code
Revision as of 17:18, 8 February 2011 by Rdm (talk | contribs) (→‎{{header|J}})
Fast 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 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

Based on j:Essays/FFT, with some simplifications, sacrificing accuracy and optimizations not visible here for clarity:

<lang j>cube =: ($~ q:@#) :. , roots =: +@]^:(0>[) (_1^2%]) ^ i.@-: rou =: ^@j.@o.@(% #)@i.@-: floop =: 4 : 'for_r. i.#$x do. (y=.{."1 y) ] x=.(+/x) ,&,:"r (-/x)*y end.' fft =: ( ] floop&.cube rou@#) :. ifft ifft =: (# %~ ] floop&.cube +@rou@#) :. fft</lang>

Example:

<lang j> require'printf'

  fmt =: [:, sprintf~&'%7.3f'@|
  (>'wave:';'fft:'),.fmt"1 (,:fft) 1 o. 2p1*3r16 * i.16

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</lang>

Note that fmt only displays the magnitude of the result, because sprintf does not support complex arguments.

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