Разбор алгоритма быстрого преобразования фурье (он же Fast Fourier transform, он же FFT).

Целью публикации является отобразить суть алгоритма, а не создание быстрой по скорости работы реализации.
Т.к. в этой программе я буду интенсивно использовать объекты вместо примитивных типов, а также рекурсию,
скорость работы будет заведомо уступать многим другим реализациям FFT, которые можно сейчас найти
в интернете в огромном количестве в открытом доступе.

В защиту данной реализации, хотел сказать что она возможно для многих будет проще для понимания 🙂

Для начала вернёмся к определению комплексных чисел:

/*
 *  Автор: Лигай Виталий vit@programmisty.com
 *  Author: VITALIY LIGAY
 */
class Complex (val x: Double, val y: Double) {

  def +(c : Complex) = new Complex(x + c.x, y + c.y)
  def -(c : Complex) = new Complex(x - c.x, y - c.y)
  def *(c : Complex) = new Complex(x * c.x - y*c.y, y*c.x + x * c.y)
  def /(c : Complex) = new Complex( (x * c.x - y*c.y)
                                   / (c.x * c.x + c.y*c.y) ,
                                   (y * c.x - x * c.y)/ (c.x * c.x + c.y*c.y))

  // sqrt(x*x + y*y)
  def abs() = math.hypot(re, im)
  def abs2() = x*x + y*y
  override def toString() = re + (if (im<0.0) "-" else "+") + math.abs(im) + "i"
}

Эта такая преамбула, теперь возьмем выжимку из русской википедии относительно принципа работы FFT.

Кстати, подробную реализацию FFT на java я нашел здесь.

/*
 *  Автор: Лигай Виталий vit@programmisty.com
 *  Author: VITALIY LIGAY
 */
import scala.math._;

object FFT {
  
  def fft(x : Array[Complex]) : Array[Complex] = {
    val N = x.length
    if (N % 2 != 0) {
      throw new RuntimeException("N должно делится на 2. В этом изюминка метода.")
    }
    // результат преобразования
    val y = new Array[Complex](N)
    if (N > 2) {
      // четная половина
      val half = new Array[Complex](N/2)
      for (n < - 0 until N/2) {
        half(n) = x(2*n)
      }
      val y0 = fft(half)
      
      // нечетная половина
      for (n <- 0 until N/2) {
        half(n)  = x(2*n + 1)
      }
      val y1 = fft(half)
      
      val ω = -2*Pi/N
      for (m <-0 until N/2) {
        // -i*2πk/N.  Особая благодарность формуле Эйлера: http://en.wikipedia.org/wiki/Euler%27s_formula
        val exp = new Complex(cos(ω*m), sin(ω*m))
        y(m)       = y0(m) + exp * y1(m)
        y(m + N/2) = y0(m) - exp * y1(m)
      }
    } else {
      // для двух элементов
      y(0) = x(0) + x(1)
      y(1) = x(0) - x(1)
    }
    y
  }

  // выводим содержимое массива 
  def show(x :Array[Complex], title: String) {
    println(title);
    println("-------------------");
    for (i <- 0 to x.length - 1) {
      printf("%f %n", x(i).abs);
    }
    println();
  }
    
  def main(a: Array[String]) = {
    val N = 16;
    val x = new Array[Complex](N);

    // я - волна, обычная волна...
    val step = 2*Pi/N;
    for (i <- 0 to N-1) {
      x(i) = new Complex(cos(i*step), sin(i*step));
    }
    show(x, "x");
    // 
    val y = fft(x);
    show(y, "y");
  }
}

Например если посчитать FFT для функции sin(t) + 0.5sin(4t), можно получить следующую картинку (красным - волна, зеленым - спектр):

P.S.: Эта первая реализация, которая пришла в голову, многое можно упростить. Будут идеи - пишите.