过河拜码头，上岸数人马。能走（苟）多远走（苟）多远，送死的人来了。
《好家伙》主题曲
从2023的国庆节到现在都龙年春节了，过去几个月，一言难尽。偶尔翻出来了上学期间写的一个小blog，看着那时候中二的文字，也是难以启齿。但也算乐于思考。往事可追，这里回顾这两个很有趣的算法，作为2024年思考的开始。
约瑟夫环问题
最早是本科数据结构课程，用链表去模拟这个问题。后来是随着Knuth教授的《具体数学》这本书的第一章用递归算式的技巧研究了这个问题，书中公式(1.10)给出了约瑟夫环问题本质是bit的循环位移。当时我看到这里的时候，立刻想到了FFT的倒序重排，两者肯定是一样的。这里不就不赘述约瑟夫问题了，强烈推荐去看Knuth教授的书，非常有趣，受益良多。
we get J(n) from n by doing a one-bit cyclic shift left! Magic
Concrete Math, Knuth
FFT算法
FFT是DFT的快速算法，DFT的公式如下:
X[k] = \sum_{n=0}^{N-1}x[n]W_{N}^{kn}, \ \ \ \ k = 0,1,\dots,N-1; W_N = e^{-j\frac{2\pi}{N}}FFT里我们认为\exists u\in \mathbf{Z}^{+} , N=2^u. 一般DSP教材里会用介绍两种办法去计算，称为时域采样和频域采样，实际上就是想办法找一个模式一样的子问题去递归。这里我们从频率抽样法入手，一步递归的本质就是找到X[2k] 和 X[2k+1] 这两个频率序列的时间序列。下面做一些简单的推理：
Y[k] = X[2k] = \sum_{n=0}^{N-1} x[n]W_{N}^{2kn} = \sum_{n=0}^{\frac{N}{2}-1}x[n]W_{N}^{2kn} + \sum_{n=\frac{N}{2}}^{N-1}x[n]W_{N}^{2kn} \\
= \sum_{n=0}^{\frac{N}{2}-1}x[n]W_{\frac{N}{2}}^{kn} + \sum_{n=0}^{\frac{N}{2}-1}x[n+\frac{N}{2}]W_{\frac{N}{2}}^{kn} = \sum_{n=0}^{\frac{N}{2}-1}(x[n]+x[n+\frac{N}{2}]) W_{\frac{N}{2}}^{kn}因此可以看出FFT基于频域采样的蝶形算法流图如下:
我们现在研究最后一排，即完成FFT后，频率下标和原始排列之间的映射关系。完全可以用Knuth的约瑟夫环问题的相同技巧。为了方便推导，我们引入一些符号：
- J(k, n): 对N=2^k个点进行FFT，上述蝶形算法后，第n个位置的频率下标, 有n\in [0,2^k)
- 比如上图，J(3, 0) = 0, J(3,1) = 4, J(3,5) = 5, J(3,6) = 3.
我们用类似Knuth的技巧，来研究J(k, n)的递归关系。共需要k轮运算。
- 当只有两个采样点的时候，k=1, 直接对应即可: J(1, n) = n, n\in[0,1]
- 对于J(k, n), n<\frac{N}{2},即上半部分，或者说，时域位置的最高级
bit是0. 这部分（位置）对应一个小规模的FFT，且频率位置是偶数，因此有 J(k, n) = 2J(k-1, n), n \in [0, \frac{N}{2}) - 对于J(k, n), n \le \frac{N}{2},即下半部分，或者说，时域位置的最高级
bit是1.这部分和上面几乎一样，就是对应奇数部分，所以有一个位置的偏移, 因此有J(k,n) = J(k, n-\frac{N}{2})+1 = J(k, n-\frac{N}{2}) + 1
利用边界条件和两个递归表达式，合计三个式子，我们就可以推演出一切：
J(1, n) = n, n\in [0, 1]\\
J(k, n) = 2J(k-1, n), n\in[0, 2^{k})\\
J(k, n) = 2J(k-1, n-2^{k})+1,n\ge 2^{k}- 上面第二个式子，n首位
bit是0，结果的末尾bit是0。 - 上面第三个式子，n首位
bit是1，结果的末尾bit是1。
利用数学归纳法，很容易可以得出，J(k,n)的结果就是把n进行bit位置倒排。
一段Erlang代码可以是：
-module(fft).
-include_lib("eunit/include/eunit.hrl").
-export([j/2]).
j(1, N) -> N;
j(K, N) when K > 1 ->
H = trunc(math:pow(2, K-1)),
case N < H of
true ->
2 * j(K-1, N);
false ->
2 * j(K-1, N - H) + 1
end.
fft_test() ->
%% k = 3, 8点fft
?assert([j(3, N) || N <- lists:seq(0, 7)] =:=
[0,4,2,6,1,5,3,7]),
%% k = 2, 4点fft
?assert([j(2, N) || N <- lists:seq(0, 3)] =:=
[0,2,1,3]),
ok.