BOB*_*BOB 5 language-agnostic algorithm bit-manipulation fft
我试图理解这个位反转算法。我找到了很多资料，但它并没有真正解释伪代码是如何工作的。例如，我从http://www.briangough.com/fftalgorithms.pdf找到了下面的伪代码
\n\nfor i = 0 ... n \xe2\x88\x92 2 do\n k = n/2\n if i < j then\n swap g(i) and g(j)\n end if\n while k \xe2\x89\xa4 j do\n j \xe2\x87\x90 j \xe2\x88\x92 k\n k \xe2\x87\x90 k/2\n end while\n j \xe2\x87\x90 j + k\nend for\n从这个伪代码来看，我不明白你为什么要这么做
\n\nswap g(i) and g(j)
当if声明是true.
另外：while循环有什么作用？如果有人可以向我解释这个伪代码，那就太好了。
下面是我在网上找到的c++代码。
\n\nvoid four1(double data[], int nn, int isign)\n{\n int n, mmax, m, j, istep, i;\n double wtemp, wr, wpr, wpi, wi, theta;\n double tempr, tempi;\n\n n = nn << 1;\n j = 1;\n for (i = 1; i < n; i += 2) {\n if (j > i) {\n tempr = data[j]; data[j] = data[i]; data[i] = tempr;\n tempr = data[j+1]; data[j+1] = data[i+1]; data[i+1] = tempr;\n }\n m = n >> 1;\n while (m >= 2 && j > m) {\n j -= m;\n m >>= 1;\n }\n j += m;\n }\n这是我发现的 FFT 源代码的完整版本
\n\n/************************************************\n* FFT code from the book Numerical Recipes in C *\n* Visit www.nr.com for the licence. *\n************************************************/\n\n// The following line must be defined before including math.h to correctly define M_PI\n#define _USE_MATH_DEFINES\n#include <math.h>\n#include <stdio.h>\n#include <stdlib.h>\n\n#define PI M_PI /* pi to machine precision, defined in math.h */\n#define TWOPI (2.0*PI)\n\n/*\n FFT/IFFT routine. (see pages 507-508 of Numerical Recipes in C)\n\n Inputs:\n data[] : array of complex* data points of size 2*NFFT+1.\n data[0] is unused,\n * the n\'th complex number x(n), for 0 <= n <= length(x)-1, is stored as:\n data[2*n+1] = real(x(n))\n data[2*n+2] = imag(x(n))\n if length(Nx) < NFFT, the remainder of the array must be padded with zeros\n\n nn : FFT order NFFT. This MUST be a power of 2 and >= length(x).\n isign: if set to 1, \n computes the forward FFT\n if set to -1, \n computes Inverse FFT - in this case the output values have\n to be manually normalized by multiplying with 1/NFFT.\n Outputs:\n data[] : The FFT or IFFT results are stored in data, overwriting the input.\n*/\n\nvoid four1(double data[], int nn, int isign)\n{\n int n, mmax, m, j, istep, i;\n double wtemp, wr, wpr, wpi, wi, theta;\n double tempr, tempi;\n\n n = nn << 1;\n j = 1;\n for (i = 1; i < n; i += 2) {\n if (j > i) {\n //swap the real part\n tempr = data[j]; data[j] = data[i]; data[i] = tempr;\n //swap the complex part\n tempr = data[j+1]; data[j+1] = data[i+1]; data[i+1] = tempr;\n }\n m = n >> 1;\n while (m >= 2 && j > m) {\n j -= m;\n m >>= 1;\n }\n j += m;\n }\n mmax = 2;\n while (n > mmax) {\n istep = 2*mmax;\n theta = TWOPI/(isign*mmax);\n wtemp = sin(0.5*theta);\n wpr = -2.0*wtemp*wtemp;\n wpi = sin(theta);\n wr = 1.0;\n wi = 0.0;\n for (m = 1; m < mmax; m += 2) {\n for (i = m; i <= n; i += istep) {\n j =i + mmax;\n tempr = wr*data[j] - wi*data[j+1];\n tempi = wr*data[j+1] + wi*data[j];\n data[j] = data[i] - tempr;\n data[j+1] = data[i+1] - tempi;\n data[i] += tempr;\n data[i+1] += tempi;\n }\n wr = (wtemp = wr)*wpr - wi*wpi + wr;\n wi = wi*wpr + wtemp*wpi + wi;\n }\n mmax = istep;\n }\n}\n\n/********************************************************\n* The following is a test routine that generates a ramp *\n* with 10 elements, finds their FFT, and then finds the *\n* original sequence using inverse FFT *\n********************************************************/\n\nint main(int argc, char * argv[])\n{\n int i;\n int Nx;\n int NFFT;\n double *x;\n double *X;\n\n /* generate a ramp with 10 numbers */\n Nx = 10;\n printf("Nx = %d\\n", Nx);\n x = (double *) malloc(Nx * sizeof(double));\n for(i=0; i<Nx; i++)\n {\n x[i] = i;\n }\n\n /* calculate NFFT as the next higher power of 2 >= Nx */\n NFFT = (int)pow(2.0, ceil(log((double)Nx)/log(2.0)));\n printf("NFFT = %d\\n", NFFT);\n\n /* allocate memory for NFFT complex numbers (note the +1) */\n X = (double *) malloc((2*NFFT+1) * sizeof(double));\n\n /* Storing x(n) in a complex array to make it work with four1. \n This is needed even though x(n) is purely real in this case. */\n for(i=0; i<Nx; i++)\n {\n X[2*i+1] = x[i];\n X[2*i+2] = 0.0;\n }\n /* pad the remainder of the array with zeros (0 + 0 j) */\n for(i=Nx; i<NFFT; i++)\n {\n X[2*i+1] = 0.0;\n X[2*i+2] = 0.0;\n }\n\n printf("\\nInput complex sequence (padded to next highest power of 2):\\n");\n for(i=0; i<NFFT; i++)\n {\n printf("x[%d] = (%.2f + j %.2f)\\n", i, X[2*i+1], X[2*i+2]);\n }\n\n /* calculate FFT */\n four1(X, NFFT, 1);\n\n printf("\\nFFT:\\n");\n for(i=0; i<NFFT; i++)\n {\n printf("X[%d] = (%.2f + j %.2f)\\n", i, X[2*i+1], X[2*i+2]);\n }\n\n /* calculate IFFT */\n four1(X, NFFT, -1);\n\n /* normalize the IFFT */\n for(i=0; i<NFFT; i++)\n {\n X[2*i+1] /= NFFT;\n X[2*i+2] /= NFFT;\n }\n\n printf("\\nComplex sequence reconstructed by IFFT:\\n");\n for(i=0; i<NFFT; i++)\n {\n printf("x[%d] = (%.2f + j %.2f)\\n", i, X[2*i+1], X[2*i+2]);\n }\n\n getchar();\n}\n\n/*\n\nNx = 10\nNFFT = 16\n\nInput complex sequence (padded to next highest power of 2):\nx[0] = (0.00 + j 0.00)\nx[1] = (1.00 + j 0.00)\nx[2] = (2.00 + j 0.00)\nx[3] = (3.00 + j 0.00)\nx[4] = (4.00 + j 0.00)\nx[5] = (5.00 + j 0.00)\nx[6] = (6.00 + j 0.00)\nx[7] = (7.00 + j 0.00)\nx[8] = (8.00 + j 0.00)\nx[9] = (9.00 + j 0.00)\nx[10] = (0.00 + j 0.00)\nx[11] = (0.00 + j 0.00)\nx[12] = (0.00 + j 0.00)\nx[13] = (0.00 + j 0.00)\nx[14] = (0.00 + j 0.00)\nx[15] = (0.00 + j 0.00)\n\nFFT:\nX[0] = (45.00 + j 0.00)\nX[1] = (-25.45 + j 16.67)\nX[2] = (10.36 + j -3.29)\nX[3] = (-9.06 + j -2.33)\nX[4] = (4.00 + j 5.00)\nX[5] = (-1.28 + j -5.64)\nX[6] = (-2.36 + j 4.71)\nX[7] = (3.80 + j -2.65)\nX[8] = (-5.00 + j 0.00)\nX[9] = (3.80 + j 2.65)\nX[10] = (-2.36 + j -4.71)\nX[11] = (-1.28 + j 5.64)\nX[12] = (4.00 + j -5.00)\nX[13] = (-9.06 + j 2.33)\nX[14] = (10.36 + j 3.29)\nX[15] = (-25.45 + j -16.67)\n\nComplex sequence reconstructed by IFFT:\nx[0] = (0.00 + j -0.00)\nx[1] = (1.00 + j -0.00)\nx[2] = (2.00 + j 0.00)\nx[3] = (3.00 + j -0.00)\nx[4] = (4.00 + j -0.00)\nx[5] = (5.00 + j 0.00)\nx[6] = (6.00 + j -0.00)\nx[7] = (7.00 + j -0.00)\nx[8] = (8.00 + j 0.00)\nx[9] = (9.00 + j 0.00)\nx[10] = (0.00 + j -0.00)\nx[11] = (0.00 + j -0.00)\nx[12] = (0.00 + j 0.00)\nx[13] = (-0.00 + j -0.00)\nx[14] = (0.00 + j 0.00)\nx[15] = (0.00 + j 0.00)\n\n*/\n位反转算法通过反转每个项目的二进制地址来创建数据集的排列；例如，在 16 项集合中，地址：
0001 0010 0011 0100 0101 0110 0111 1000 1001 1010 1011 1100 1101 1110 1111
将更改为：
1000 0100 1100 0010 1010 0110 1110 0001 1001 0101 1101 0011 1011 0111 1111
然后相应的项将移动到新地址。
或者用十进制表示：
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
变为
0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15
伪代码中的 while 循环的作用是将变量 j 设置为该序列。（顺便说一句，j 的初始值应该为 0）。
您将看到该序列的组成如下：
0
0 1
0 2 1 3
0 4 2 6 1 5 3 7
0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15
每个序列都是通过将前一个版本乘以 2，然后添加 1 来重复得到的。或者以另一种方式看待它：通过重复前面的序列，与值 + n/2 交错（这更接近地描述了算法中发生的情况）。
0
0 1
0 2 1 3
0 4 2 6 1 5 3 7
0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15
然后，i 和 j 项会在 for 循环的每次迭代中交换，但前提是 i < j；否则，每个项目都将被交换到新位置（例如，当 i = 3 且 j = 12 时），然后再次返回（当 i = 12 且 j = 3 时）。
0
0 1
0 2 1 3
0 4 2 6 1 5 3 7
0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15
您发现的 C++ 代码似乎使用序列的对称性以双步循环遍历它。但它不会产生正确的结果，所以要么是一次失败的尝试，要么它的设计目的是完全不同的。这是使用两步思想的版本：
function bitReversal(data) {
var n = data.length;
var j = 0;
for (i = 0; i < n - 1; i++) {
var k = n / 2;
if (i < j) {
var temp = data[i]; data[i] = data[j]; data[j] = temp;
}
while (k <= j) {
j -= k;
k /= 2;
}
j += k;
}
return(data);
}
console.log(bitReversal([0,1]));
console.log(bitReversal([0,1,2,3]));
console.log(bitReversal([0,1,2,3,4,5,6,7]));
console.log(bitReversal([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]));
console.log(bitReversal(["a","b","c","d","e","f","g","h","i","j","k","l","m","n","o","p"]));