Luiz-Monad · January 18, 2021 00:28
diff --git a/fft.apl b/fft.apl
 [Math/APL] (warning: walls of text and equations)

 I loved APL !
 Exponential negativo em alpha (primeiro parametro) te dá numeros complexos ! ( http://www.jsoftware.com/papers/satn40.htm )
 Geralmente eu odeio trabalhar com arrays em C, é muito cumbersome. Mas em APL arrays são lindos.
 Mas o mais legal mesmo é a implementação das funções trigonometricas diadicas.
 This language is so good for simple mathematics, now I know why such a beaultiful language didn't work in the boring 'real life'.

 Ontem estava estudando a transformada de Fourier, dai falaram de APL, dai eu tive um sonho.
 Hum, lets implement it !

 ⍝ Decimation-in-time radix-2 ( Cooley–Tukey )
 ⍝ EqI: 
 ⍝ $X_k = \sum_{n=0}^{N-1} x_n e^{ \frac{ -2 \pi i }{N} nk }$ (fig0)

 ⍝ Eu não vou colocar todos os passos intermediarios, assista esse video ( https://www.youtube.com/watch?v=htCj9exbGo0 ), I ignored the 'C' bullshit because it didn't help, also ( https://www.youtube.com/watch?v=1mVbZLHLaf0 ) helps too.

 ⍝ Separar entre pares e impares.
 ⍝ $X_k = \sum_{m=0}^{N/2-1} x_{2m} e^{ \frac{ -2 \pi i }{N} (2m)k } + \sum_{m=0}^{N/2-1} x_{2m+1} e^{ \frac{ -2 \pi i }{N} (2m+1)k }$ (fig1)

 ⍝ Extrair fator comum.
 ⍝ EqII:
 ⍝ $e^{ \frac{ -2 \pi i }{N} k }$ (fig2)

 ⍝ Mas antes vamos converter o fator comum para APL. (APL se le da direita para esquerda e não tem precedencia de operador, exceto parens)  
 ⍝ e^{ \frac{ -2 \pi i }{N} k }
      * (⍟¯1) × ¯2 × k÷N
 ⍝ Lets test with k=1 and N=32
 0.980785J¯0.19509 
 ⍝ Exactly what Wolfram tell us ( https://www.wolframalpha.com/input/?i=e%5E(+(-2+%5Cpi+i)+%2F(32)) )

 ⍝ ⍟¯1 é o mesmo que −2πi. 
 ⍝ -2π é ¯2×○1. 
 ⍝ ⍟ é o logaritimo natural.
 ⍝ ○ é a função circular trigonometrica. ○1 é π.
 ⍝ Oh, eu nunca tinha percebido isto antes na Matematica, makes totally sense! (proof left as an exercise for the reader)

 ⍝ What we got so far.
 ⍝ $X_k = \sum_{m=0}^{N/2-1} x_{2m} e^{ \frac{ -2 \pi i }{N} (2m)k } + \sum_{m=0}^{N/2-1} x_{2m+1} e^{ \frac{ -2 \pi i }{N} (2m+1)k }$ (fig3)
 ⍝ Extracting commom factor.
 ⍝ EqIII:
 ⍝ $X_k = \sum_{m=0}^{N/2-1} x_{2m} e^{ \frac{ -2 \pi i }{N/2} mk } + e^{ \frac{ -2 \pi i }{N} k } \sum_{m=0}^{N/2-1} x_{2m+1} e^{ \frac{ -2 \pi i }{N/2} mk }$ (fig4)
 ⍝ Lets simplify (replacing from Eq I).
 ⍝ $X_k = Even_k + e^{ \frac{ -2 \pi i }{N} k } Odd_k $ (fig5)
 ⍝ Even_k and Odd_k are recursive FFTs themselves (recursion and FP, take that, divide and conquer!)
 ⍝ Now we only have simple equations, and array accesses, nice!

 ⍝ As per niquist theorem we know that we only need half of the periodicity as the function will repeat values.
 ⍝ $Even_{k+\frac{N}{2}} = Even_k$ (fig6) (Odd is exactly the same)
 ⍝ Então o N pode ser reduzido para N/2.

 ⍝ As equações ficam da seguinte forma.
 ⍝ EqIV:
 ⍝ $\newline X_k = Even_k + factor * Odd_k \: for \, 1 <= k <= N/2 
 ⍝  \newline X_k = Even_{k-N/2} + factor * Odd_{(k-N/2)} \: for \, N/2 <= k <= N $ (fig7)
 ⍝ $Even_k$ e $Odd_k$ são recursive calls em EqI, pois foi substituidos em EqIII.

 ⍝ Agora vamos usar a simetria ( https://youtu.be/1mVbZLHLaf0?t=516 )
 ⍝ Então plugando $k+\frac{N}{2}$ na EqII e sabendo da simetria (tiddle factor).
 ⍝ $e^{ \frac{ -2 \pi i }{N} (k+\frac{N}{2}) }$ (fig8)
 ⍝ A equação fica:
 ⍝ EqV:
 ⍝ $-e^{ \frac{ -2 \pi i }{N} k }$ (fig9)
 ⍝ Substituindo devolta EqV em EqIV.
 ⍝ $\newline X_k = Even_k + factor * Odd_k
 ⍝  \newline X_k = Even_k - factor * Odd_k$ (fig10)

 ⍝ Agora podemos enfiar em APL:
 ⍝ the factor
      f ← * (⍟¯1) × ¯2 × k÷N
 ⍝ the recursion
      T ← FFT k-N÷2
 ⍝ equation
 ⍝      O[k] ← T[k] + f × T[k+N÷2]      ⍝ even
 ⍝      O[k+N÷2] ← T[k] - f × T[k+N÷2]  ⍝ odd

 ⍝ Lets test it.
      X ← ⍳8
      s ← 8
      N ← 8
      k ← 1
      f ← * (⍟¯1) × ¯2 × k÷N
      T[k] + f × T[k+N÷2]
 4.53553J¯3.53553
      T[k] - f × T[k+N÷2]
 ¯2.53553J3.53553
 ⍝ Perfect! ( https://www.wolframalpha.com/input/?i=1%2B(e%5E(+(+-2+%5Cpi+i+)%2F(8)+1+))*5 )

 Agora vamos escrever para usar os primitivos de array de APL.

 ⍝ Primeiro precisamos cortar a sequencia em impar e par e depois aplicar a function recursivamente.
 ⍝ Reshape X to be N÷2 rows per 2 cols matrix.
      T ← X ⍴⍨ 2 ,⍨ N÷2

 ⍝ then apply columns into FFT.
      FFT ← {⍵} ⍝ dummy
      Even ← FFT 1⌷[2]T  ⍝ or T[;1] 
      Odd ← FFT 2⌷[2]T   ⍝ or T[;2] 
 ⍝ combined
      EvenOdd ← Even Odd

 ⍝ Dai executa a somatoria de k.
 ⍝ ⍺ is even, ⍵ is odd
 ⍝ ⍺ ← T[k] ⋄ ⍵ ← T[k+N÷2]
 ⍝ for every row do thing inside {}.
 ⍝ we just replaced the variables
      Y1 ← {⍺ + f × ⍵} / EvenOdd
      Y2 ← {⍺ - f × ⍵} / EvenOdd
 ⍝ The final result is Y1 Y2 combined.
      Y ← Y1 Y2

 ⍝ O resultado são 4 numeros complexos para Y1 e mais 4 para Y2.
 ⍝ ┌───────────────────────────────────────────────────────────────────┐
 ⍝ │2.41421J¯1.41421 5.82843J¯2.82843 9.24264J¯4.24264 12.6569J¯5.65685│
 ⍝ └───────────────────────────────────────────────────────────────────┘
 ⍝ ┌───────────────────────────────────────────────────────────────────┐
 ⍝ │¯0.414214J1.41421 0.171573J2.82843 0.757359J4.24264 1.34315J5.65685│
 ⍝ └───────────────────────────────────────────────────────────────────┘

 ⍝ Vamos fazer um curto circuito agora.
 ⍝ apply reduce to union then and map FFT with dieresis.
      EvenOdd ← FFT ¨ ∪ ⌿ T

 ⍝ Lets fix a mistake, f deveria ser o termo da EqII e o k é parametro.
 ⍝ We also need the index given
      omg ← { * (⍟¯1) × ¯2 × ⍵÷N }
      EvenOdd ← FFT ¨ ∪ ⌿ T
      Ix ← { ⍵ (⍳⍴⍵ ⍴⍨ ⍴⍵) } ¨ EvenOdd
      Y1 ← { ⍺ + ⍵ × omg ⍵ } / {⍵ Ix} ¨ EvenOdd
      Y2 ← { ⍺ - ⍵ × omg ⍵ } / {⍵ Ix} ¨ EvenOdd

 ⍝ So far
      X ← ⍳8
      s ← 8
      N ← 8
      k ← 1
      omg ← { * (⍟¯1) × ¯2 × ⍵÷N }
      T ← X ⍴⍨ 2 ,⍨ N÷2
      FFT ← {⍵} 
      EvenOdd ← FFT ¨ ∪ ⌿ T
      Y1 ← { ⍺ + ⍵ × omg k } / EvenOdd
      Y2 ← { ⍺ - ⍵ × omg k } / EvenOdd

 ⍝ Agora só construir a recursao.

 FFT ← { ⍺ = 1 : ⍵
  f ← * (⍟¯1) × ¯2 × k÷N
  T ← X ⍴⍨ 2 ,⍨ N÷2
  EvenOdd ← ∇ ¨ ∪ ⌿ T
  Y1 ← {⍺ + f × ⍵} / EvenOdd
  Y2 ← {⍺ - f × ⍵} / EvenOdd

 (⍺-1)∇ ⍺⍺ ⍵ 
 }

 ⍝ Todo: butterfly optimization.
 ⍝ Todo: using matrices instead of iteration ( https://youtu.be/1mVbZLHLaf0?t=3348 ) .
 ⍝ Todo: more optimization after reading this ( http://cnx.org/contents/ulXtQbN7@15/Implementing-FFTs-in-Practice )

 Tudo porque eu estava estudando Quantum Physics e precisei de reestudar Calculo, este video aqui foi a inspiração ( https://www.youtube.com/watch?v=r18Gi8lSkfM ). (There's crazy people like me studying Calculus for fun) (Why study only one thing when you can study 2 or 3 things at the same time)

 Proximo episodio: Eigenvalues ! just kidding
	[Math/APL] (warning: walls of text and equations)

	I loved APL !
	Exponential negativo em alpha (primeiro parametro) te dá numeros complexos ! ( http://www.jsoftware.com/papers/satn40.htm )
	Geralmente eu odeio trabalhar com arrays em C, é muito cumbersome. Mas em APL arrays são lindos.
	Mas o mais legal mesmo é a implementação das funções trigonometricas diadicas.
	This language is so good for simple mathematics, now I know why such a beaultiful language didn't work in the boring 'real life'.

	Ontem estava estudando a transformada de Fourier, dai falaram de APL, dai eu tive um sonho.
	Hum, lets implement it !

	⍝ Decimation-in-time radix-2 ( Cooley–Tukey )
	⍝ EqI:
	⍝ $X_k = \sum_{n=0}^{N-1} x_n e^{ \frac{ -2 \pi i }{N} nk }$ (fig0)

	⍝ Eu não vou colocar todos os passos intermediarios, assista esse video ( https://www.youtube.com/watch?v=htCj9exbGo0 ), I ignored the 'C' bullshit because it didn't help, also ( https://www.youtube.com/watch?v=1mVbZLHLaf0 ) helps too.

	⍝ Separar entre pares e impares.
	⍝ $X_k = \sum_{m=0}^{N/2-1} x_{2m} e^{ \frac{ -2 \pi i }{N} (2m)k } + \sum_{m=0}^{N/2-1} x_{2m+1} e^{ \frac{ -2 \pi i }{N} (2m+1)k }$ (fig1)

	⍝ Extrair fator comum.
	⍝ EqII:
	⍝ $e^{ \frac{ -2 \pi i }{N} k }$ (fig2)

	⍝ Mas antes vamos converter o fator comum para APL. (APL se le da direita para esquerda e não tem precedencia de operador, exceto parens)
	⍝ e^{ \frac{ -2 \pi i }{N} k }
	* (⍟¯1) × ¯2 × k÷N
	⍝ Lets test with k=1 and N=32
	0.980785J¯0.19509
	⍝ Exactly what Wolfram tell us ( https://www.wolframalpha.com/input/?i=e%5E(+(-2+%5Cpi+i)+%2F(32)) )

	⍝ ⍟¯1 é o mesmo que −2πi.
	⍝ -2π é ¯2×○1.
	⍝ ⍟ é o logaritimo natural.
	⍝ ○ é a função circular trigonometrica. ○1 é π.
	⍝ Oh, eu nunca tinha percebido isto antes na Matematica, makes totally sense! (proof left as an exercise for the reader)

	⍝ What we got so far.
	⍝ $X_k = \sum_{m=0}^{N/2-1} x_{2m} e^{ \frac{ -2 \pi i }{N} (2m)k } + \sum_{m=0}^{N/2-1} x_{2m+1} e^{ \frac{ -2 \pi i }{N} (2m+1)k }$ (fig3)
	⍝ Extracting commom factor.
	⍝ EqIII:
	⍝ $X_k = \sum_{m=0}^{N/2-1} x_{2m} e^{ \frac{ -2 \pi i }{N/2} mk } + e^{ \frac{ -2 \pi i }{N} k } \sum_{m=0}^{N/2-1} x_{2m+1} e^{ \frac{ -2 \pi i }{N/2} mk }$ (fig4)
	⍝ Lets simplify (replacing from Eq I).
	⍝ $X_k = Even_k + e^{ \frac{ -2 \pi i }{N} k } Odd_k $ (fig5)
	⍝ Even_k and Odd_k are recursive FFTs themselves (recursion and FP, take that, divide and conquer!)
	⍝ Now we only have simple equations, and array accesses, nice!

	⍝ As per niquist theorem we know that we only need half of the periodicity as the function will repeat values.
	⍝ $Even_{k+\frac{N}{2}} = Even_k$ (fig6) (Odd is exactly the same)
	⍝ Então o N pode ser reduzido para N/2.

	⍝ As equações ficam da seguinte forma.
	⍝ EqIV:
	⍝ $\newline X_k = Even_k + factor * Odd_k \: for \, 1 <= k <= N/2
	⍝ \newline X_k = Even_{k-N/2} + factor * Odd_{(k-N/2)} \: for \, N/2 <= k <= N $ (fig7)
	⍝ $Even_k$ e $Odd_k$ são recursive calls em EqI, pois foi substituidos em EqIII.

	⍝ Agora vamos usar a simetria ( https://youtu.be/1mVbZLHLaf0?t=516 )
	⍝ Então plugando $k+\frac{N}{2}$ na EqII e sabendo da simetria (tiddle factor).
	⍝ $e^{ \frac{ -2 \pi i }{N} (k+\frac{N}{2}) }$ (fig8)
	⍝ A equação fica:
	⍝ EqV:
	⍝ $-e^{ \frac{ -2 \pi i }{N} k }$ (fig9)
	⍝ Substituindo devolta EqV em EqIV.
	⍝ $\newline X_k = Even_k + factor * Odd_k
	⍝ \newline X_k = Even_k - factor * Odd_k$ (fig10)

	⍝ Agora podemos enfiar em APL:
	⍝ the factor
	f ← * (⍟¯1) × ¯2 × k÷N
	⍝ the recursion
	T ← FFT k-N÷2
	⍝ equation
	⍝ O[k] ← T[k] + f × T[k+N÷2] ⍝ even
	⍝ O[k+N÷2] ← T[k] - f × T[k+N÷2] ⍝ odd

	⍝ Lets test it.
	X ← ⍳8
	s ← 8
	N ← 8
	k ← 1
	f ← * (⍟¯1) × ¯2 × k÷N
	T[k] + f × T[k+N÷2]
	4.53553J¯3.53553
	T[k] - f × T[k+N÷2]
	¯2.53553J3.53553
	⍝ Perfect! ( https://www.wolframalpha.com/input/?i=1%2B(e%5E(+(+-2+%5Cpi+i+)%2F(8)+1+))*5 )

	Agora vamos escrever para usar os primitivos de array de APL.

	⍝ Primeiro precisamos cortar a sequencia em impar e par e depois aplicar a function recursivamente.
	⍝ Reshape X to be N÷2 rows per 2 cols matrix.
	T ← X ⍴⍨ 2 ,⍨ N÷2

	⍝ then apply columns into FFT.
	FFT ← {⍵} ⍝ dummy
	Even ← FFT 1⌷[2]T ⍝ or T[;1]
	Odd ← FFT 2⌷[2]T ⍝ or T[;2]
	⍝ combined
	EvenOdd ← Even Odd

	⍝ Dai executa a somatoria de k.
	⍝ ⍺ is even, ⍵ is odd
	⍝ ⍺ ← T[k] ⋄ ⍵ ← T[k+N÷2]
	⍝ for every row do thing inside {}.
	⍝ we just replaced the variables
	Y1 ← {⍺ + f × ⍵} / EvenOdd
	Y2 ← {⍺ - f × ⍵} / EvenOdd
	⍝ The final result is Y1 Y2 combined.
	Y ← Y1 Y2

	⍝ O resultado são 4 numeros complexos para Y1 e mais 4 para Y2.
	⍝ ┌───────────────────────────────────────────────────────────────────┐
	⍝ │2.41421J¯1.41421 5.82843J¯2.82843 9.24264J¯4.24264 12.6569J¯5.65685│
	⍝ └───────────────────────────────────────────────────────────────────┘
	⍝ ┌───────────────────────────────────────────────────────────────────┐
	⍝ │¯0.414214J1.41421 0.171573J2.82843 0.757359J4.24264 1.34315J5.65685│
	⍝ └───────────────────────────────────────────────────────────────────┘

	⍝ Vamos fazer um curto circuito agora.
	⍝ apply reduce to union then and map FFT with dieresis.
	EvenOdd ← FFT ¨ ∪ ⌿ T

	⍝ Lets fix a mistake, f deveria ser o termo da EqII e o k é parametro.
	⍝ We also need the index given
	omg ← { * (⍟¯1) × ¯2 × ⍵÷N }
	EvenOdd ← FFT ¨ ∪ ⌿ T
	Ix ← { ⍵ (⍳⍴⍵ ⍴⍨ ⍴⍵) } ¨ EvenOdd
	Y1 ← { ⍺ + ⍵ × omg ⍵ } / {⍵ Ix} ¨ EvenOdd
	Y2 ← { ⍺ - ⍵ × omg ⍵ } / {⍵ Ix} ¨ EvenOdd

	⍝ So far
	X ← ⍳8
	s ← 8
	N ← 8
	k ← 1
	omg ← { * (⍟¯1) × ¯2 × ⍵÷N }
	T ← X ⍴⍨ 2 ,⍨ N÷2
	FFT ← {⍵}
	EvenOdd ← FFT ¨ ∪ ⌿ T
	Y1 ← { ⍺ + ⍵ × omg k } / EvenOdd
	Y2 ← { ⍺ - ⍵ × omg k } / EvenOdd

	⍝ Agora só construir a recursao.

	FFT ← { ⍺ = 1 : ⍵
	f ← * (⍟¯1) × ¯2 × k÷N
	T ← X ⍴⍨ 2 ,⍨ N÷2
	EvenOdd ← ∇ ¨ ∪ ⌿ T
	Y1 ← {⍺ + f × ⍵} / EvenOdd
	Y2 ← {⍺ - f × ⍵} / EvenOdd

	(⍺-1)∇ ⍺⍺ ⍵
	}

	⍝ Todo: butterfly optimization.
	⍝ Todo: using matrices instead of iteration ( https://youtu.be/1mVbZLHLaf0?t=3348 ) .
	⍝ Todo: more optimization after reading this ( http://cnx.org/contents/ulXtQbN7@15/Implementing-FFTs-in-Practice )

	Tudo porque eu estava estudando Quantum Physics e precisei de reestudar Calculo, este video aqui foi a inspiração ( https://www.youtube.com/watch?v=r18Gi8lSkfM ). (There's crazy people like me studying Calculus for fun) (Why study only one thing when you can study 2 or 3 things at the same time)

	Proximo episodio: Eigenvalues ! just kidding