我正在尝试确定一种“简单”的方法来计算 FFT 的哪些输入需要在其各个阶段一起“蝴蝶”。我正在看这样的图表：
在点 FFT 图中，
我希望对于 FFT 的每个阶段，我可以有一个从到（FFT 长度）的简单计数器，并且我可以进行非常便宜的位操作，以正确的顺序形成正确的索引。
第 1 阶段的模式相当明显：有一个从 0 到 15 的 4 位计数器，并完全反转计数器中的位：
我相信第 2 阶段是这样的：位反转计数器，然后交换最左边的 2 位。
是否有适用于任何 FFT 任何阶段的通用公式？
假设。所以计数器选择到之间的值，并且有个阶段。根据您的评论，阶段的输入-输出关系可以由下式给出：
位（此处为 4 位）表示二进制形式的计数器。
将计数器向右循环移动
位反转结果。
以十进制表示。
下面的 Matlab 代码可以给出向量形式的输出，如注释中所述：
N = 16;
stage = 1 %Do it here for stage #1
c = 0:N-1;
c_bin = de2bi(c);
butterfly_idx = bi2de(fliplr(circshift(c_bin',stage-1)'))
N=16和 stage=2`的示例：
c_bin =
0 0 0 0
1 0 0 0
0 1 0 0
1 1 0 0
0 0 1 0
1 0 1 0
0 1 1 0
1 1 1 0
0 0 0 1
1 0 0 1
0 1 0 1
1 1 0 1
0 0 1 1
1 0 1 1
0 1 1 1
1 1 1 1
fliplr(circshift(c_bin',1)')
ans =
0 0 0 0
0 0 1 0
0 1 0 0
0 1 1 0
1 0 0 0
1 0 1 0
1 1 0 0
1 1 1 0
0 0 0 1
0 0 1 1
0 1 0 1
0 1 1 1
1 0 0 1
1 0 1 1
1 1 0 1
1 1 1 1
bi2de(fliplr(circshift(c_bin',1)'))
ans =
0
4
2
6
1
5
3
7
8
12
10
14
9
13
11
15
这是第二阶段的预期结果。
因此，这个来自 loooong 之前的 FFT 代码非常简单，它实现了频率抽取 (FFT) 和时间抽取 (iFFT)。我认为您将能够找到与您的 N=16 FFT 图匹配的一段代码。
/*
A set of utility programs to compute the Fast Fourier Transform (FFT):
N-1
X[k] = SUM { x[n]exp(-j2 pi nk/N) }
n=0
and inverse Fast Fourier Transform (iFFT):
N-1
x[n] = 1/N SUM { X[k]exp(+j2 pi nk/N) }
k=0
To speed things up, multiplication by 1 and j are avoided. The input, x[], is an array of
complex numbers (pairs of doubles) of length N = 2^n. The calling program supplies
n = log2(N) not the array length, N. The output is placed in BIT REVERSED order in x[].
Call bit_reverse(x, n) to swap back to normal order, if needed. However, iFFT(X, n, stbl)
requires its input, X[], to be in bit reversed order making bit reversing unnecessary in
some cases, such as convolution. stbl is an array of doubles of length N/4 containing the
sine function from 0 to pi/2 used to compute the FFT. Call sin_table(stbl, n) ONLY ONCE
before either FFT(x, n, stbl) or iFFT(X, n, stbl) to set up a sin table for FFT computation.
Written ca. 1985 in THINK C by Robert Bristow-Johnson.
*/
#define HALFPI 1.570796326794897
#define PI 3.141592653589793
#define TWOPI 6.283185307179586
// #include "complex.h"
typedef struct {
double real;
double imag;
} complex;
#define Re(z) (z).real
#define Im(z) (z).imag
void FFT(complex *x, int n, double *stbl)
{
long size;
register long length, step, stepsize, end;
register complex *top, *bottom; /* top & bottom of FFT butterfly */
complex temp;
size = 1L<<(n-2);
end = (long)x + 4*sizeof(temp)*size;
length = size;
stepsize = 1;
while ( length >= 1)
{
top = x;
while ((long)top < end)
{
bottom = top + 2*length;
Re(temp) = Re(*top) - Re(*bottom); /* butterfly: twiddle = 1 */
Im(temp) = Im(*top) - Im(*bottom);
Re(*top) += Re(*bottom);
Im(*top) += Im(*bottom);
*bottom = temp;
top++;
bottom++;
for (step = stepsize; step < size; step += stepsize)
{
Re(temp) = Re(*top) - Re(*bottom); /* butterfly: twiddle = exp(-j theta) */
Im(temp) = Im(*top) - Im(*bottom);
Re(*top) += Re(*bottom);
Im(*top) += Im(*bottom);
Re(*bottom) = Re(temp)*stbl[size - step] + Im(temp)*stbl[step];
Im(*bottom) = Im(temp)*stbl[size - step] - Re(temp)*stbl[step];
top++;
bottom++;
}
Re(temp) = Im(*top) - Im(*bottom); /* butterfly: twiddle = -j */
Im(temp) = Re(*bottom) - Re(*top);
Re(*top) += Re(*bottom);
Im(*top) += Im(*bottom);
*bottom = temp;
top++;
bottom++;
for (step = stepsize; step < size; step += stepsize)
{
Re(temp) = Im(*top) - Im(*bottom); /* butterfly: twiddle = -j*exp(-j theta) */
Im(temp) = Re(*bottom) - Re(*top);
Re(*top) += Re(*bottom);
Im(*top) += Im(*bottom);
Re(*bottom) = Re(temp)*stbl[size - step] + Im(temp)*stbl[step];
Im(*bottom) = Im(temp)*stbl[size - step] - Re(temp)*stbl[step];
top++;
bottom++;
}
top = bottom;
}
length >>= 1;
stepsize <<= 1;
}
top = x;
bottom = x + 1;
while ((long)top < end)
{
Re(temp) = Re(*top) - Re(*bottom); /* butterfly: twiddle = 1 */
Im(temp) = Im(*top) - Im(*bottom);
Re(*top) += Re(*bottom);
Im(*top) += Im(*bottom);
*bottom = temp;
top += 2;
bottom += 2;
}
}
void iFFT(complex *X, int n, double *stbl)
{
long size;
register long length, step, stepsize, end;
double scale;
register complex *top, *bottom; /* top & bottom of FFT butterfly */
complex temp;
size = 1L<<(n-2);
end = (long)X + 4*sizeof(temp)*size;
scale = 0.25/size;
top = X;
bottom = X + 1;
while ((long)top < end)
{
Re(temp) = (Re(*top) - Re(*bottom))*scale; /* butterfly: twiddle = 1/N */
Im(temp) = (Im(*top) - Im(*bottom))*scale;
Re(*top) = (Re(*top) + Re(*bottom))*scale;
Im(*top) = (Im(*top) + Im(*bottom))*scale;
*bottom = temp;
top += 2;
bottom += 2;
}
length = 1;
stepsize = size;
while ( stepsize >= 1)
{
top = X;
while ((long)top < end)
{
bottom = top + 2*length;
temp = *bottom; /* butterfly: twiddle = 1 */
Re(*bottom) = Re(*top) - Re(temp);
Im(*bottom) = Im(*top) - Im(temp);
Re(*top) += Re(temp);
Im(*top) += Im(temp);
top++;
bottom++;
for (step = stepsize; step < size; step += stepsize)
{ /* butterfly: twiddle = exp(+j theta) */
Re(temp) = Re(*bottom)*stbl[size - step] - Im(*bottom)*stbl[step];
Im(temp) = Im(*bottom)*stbl[size - step] + Re(*bottom)*stbl[step];
Re(*bottom) = Re(*top) - Re(temp);
Im(*bottom) = Im(*top) - Im(temp);
Re(*top) += Re(temp);
Im(*top) += Im(temp);
top++;
bottom++;
}
Re(temp) = -Im(*bottom); /* butterfly: twiddle = +j */
Im(temp) = Re(*bottom);
Re(*bottom) = Re(*top) - Re(temp);
Im(*bottom) = Im(*top) - Im(temp);
Re(*top) += Re(temp);
Im(*top) += Im(temp);
top++;
bottom++;
for (step = stepsize; step < size; step += stepsize)
{ /* butterfly: twiddle = +j*exp(+j theta) */
Re(temp) = -Im(*bottom)*stbl[size - step] - Re(*bottom)*stbl[step];
Im(temp) = Re(*bottom)*stbl[size - step] - Im(*bottom)*stbl[step];
Re(*bottom) = Re(*top) - Re(temp);
Im(*bottom) = Im(*top) - Im(temp);
Re(*top) += Re(temp);
Im(*top) += Im(temp);
top++;
bottom++;
}
top = bottom;
}
length <<= 1;
stepsize >>= 1;
}
}
void sin_table(double *stbl, int n)
{
register long size, i;
double theta;
size = 1L<<(n-2);
theta = HALFPI/size;
for (i = 0; i < size; i++)
{
stbl[i] = sin(theta*i);
}
}
void bit_reverse(register complex *x, int n)
{
complex temp;
register long k, i, r, size, count;
size = (1L<<n) - 1L;
for (k = 1L; k < size; k++)
{
i = k;
r = 0;
for (count = n; count > 0; count--)
{
r <<= 1;
r += (i & 0x00000001L);
i >>= 1;
}
if (r > k)
{
temp = x[r];
x[r] = x[k];
x[k] = temp;
}
}
}