我的代码是矩阵的非常直接的实现，公式，所以我不清楚问题出在什么地方
// |
// radix-4 butterfly matrix form | complex multiplication
// |
// +- -+ +- -+ | a+ib
// X[0] = | 1 1 1 1 | |x[0]| | * c+id
// X[1] = | 1 -i -1 i | |x[1]| | -------
// X[2] = | 1 -1 1 -1 | |x[2]| | ac + ibc
// X[3] = | 1 i -1 -i | |x[3]| | iad - bd
// +- -+ +- -+ | ------------------
// | (ac-bd) + i(bc+ad)
// |
有人可以发现我出了错吗？
，谢谢，
-David
typedef double fp; // base floating-point type
// naiive N-point DFT implementation as reference to check fft implementation against
//
void dft(int inv, struct cfp *x, struct cfp *y, int N) {
long int i, j;
struct cfp w;
fp ang;
for(i=0; i<N; i++) { // do N-point FFT/IFFT
y[i].r = y[i].i = 0;
if (inv) ang = 2*PI*(fp)i/(fp)N;
else ang = -2*PI*(fp)i/(fp)N;
for (j=0; j<N; j++) {
w.r = cos(j*ang);
w.i = sin(j*ang);
y[i].r += (x[j].r * w.r - x[j].i * w.i);
y[i].i += (x[j].r * w.i + x[j].i * w.r);
}
}
// scale output in the case of an IFFT
if (inv) {
for (i=0; i<N; i++) {
y[i].r = y[i].r/(fp)N;
y[i].i = y[i].i/(fp)N;
}
}
} // dft()
void r4fft4(int inv, int reorder, struct cfp *x, struct cfp *y) {
struct cfp x1[4], w[4];
fp ang, temp;
int i;
// |
// radix-4 butterfly matrix form | complex multiplication
// |
// +- -+ +- -+ | a+ib
// y[0] = | 1 1 1 1 | |x[0]| | * c+id
// y[1] = | 1 -i -1 i | |x[1]| | -------
// y[2] = | 1 -1 1 -1 | |x[2]| | ac + ibc
// y[3] = | 1 i -1 -i | |x[3]| | iad - bd
// +- -+ +- -+ | ------------------
// | (ac-bd) + i(bc+ad)
// |
if (inv) ang = 2*PI/(fp)4; // invert sign for IFFT
else ang = -2*PI/(fp)4;
//
w[1].r = cos(ang*1); w[1].i = sin(ang*1); // twiddle1 = exp(-2*pi/4 * 1);
w[2].r = cos(ang*2); w[2].i = sin(ang*2); // twiddle2 = exp(-2*pi/4 * 2);
w[3].r = cos(ang*3); w[3].i = sin(ang*3); // twiddle3 = exp(-2*pi/4 * 3);
// *1 *1 *1 *1
y[0].r = x[0].r + x[1].r + x[2].r + x[3].r;
y[0].i = x[0].i + x[1].i + x[2].i + x[3].i;
// *1 *-i *-1 *i
x1[1].r = x[0].r + x[1].i - x[2].r - x[3].i;
x1[1].i = x[0].i - x[1].r - x[2].i + x[3].r;
// *1 *-1 *1 *-1
x1[2].r = x[0].r - x[1].r + x[2].r - x[3].r;
x1[2].i = x[0].i - x[1].i + x[2].i - x[3].i;
// *1 *i *-1 *-i
x1[3].r = x[0].r - x[1].i - x[2].r + x[3].i;
x1[3].i = x[0].i + x[1].r - x[2].i - x[3].r;
//
y[1].r = x1[1].r*w[1].r - x1[1].i*w[1].i; // scale radix-4 output
y[1].i = x1[1].i*w[1].r + x1[1].r*w[1].i;
//
y[2].r = x1[2].r*w[2].r - x1[2].i*w[2].i; // scale radix-4 output
y[2].i = x1[2].i*w[2].r + x1[2].r*w[2].i;
//
y[3].r = x1[3].r*w[3].r - x1[3].i*w[3].i; // scale radix-4 output
y[3].i = x1[3].i*w[3].r + x1[3].r*w[3].i;
// reorder output stage ... mystery as to why I need this
if (reorder) {
temp = y[1].r;
y[1].r = -1*y[1].i;
y[1].i = temp;
//
y[2].r = -1*y[2].r;
//
temp = y[3].r;
y[3].r = y[3].i;
y[3].i = -1*temp;
}
// scale output for inverse FFT
if (inv) {
for (i=0; i<4; i++) { // scale output by 1/N for IFFT
y[i].r = y[i].r/(fp)4;
y[i].i = y[i].i/(fp)4;
}
}
} // r4fft4()
#1 楼
我刚刚将S. Burrus Fortran代码的radix-4 DIF fft移植到了Java。实际上，它缺乏一些优化，首先是表驱动的旋转因子（正弦和余弦因子应预先计算）。这样可以使fft更快一些（也许50％）。为此，我必须有所作为，但是如果有人给出正确答案，我将非常高兴和感激。我将尽快发布优化的代码，希望可以通过一些速度测试与radix-2算法进行比较。
此外，不会删除与1和sqrt（-1）的乘法。删除它们会加快速度。
但是整体IMHO基数4似乎比基数2快不超过25％，所以我不知道速度/复杂度比率是否真的值得。
请记住，像FFTW这样经过非常优化的库在很大程度上已经得到使用，因此，这可能只是个人的“多元化”！将其移植到C，C ++或C＃应该非常容易。
public static void FFTR4(double[] X, double[] Y, int N, int M) {
// N = 4 ^ M
int N1,N2;
int I1, I2, I3;
double CO1,CO2,CO3,SI1,SI2,SI3;
double A,B,C,E;
double R1,R2,R3,R4;
double S1,S2,S3,S4;
// N = 1 << (M+M);
N2 = N;
I2 = 0; I3 = 0;
for (int K=0; K<M; ++K) {
N1 = N2;
N2 = N2 / 4;
E = PI2 / (double)N1;
A = 0.0;
for (int J=0; J < N2; ++J) {
A = J*E;
B = A + A;
C = A + B;
//Should be pre-calculated for optimization
CO1 = Math.cos(A);
CO2 = Math.cos(B);
CO3 = Math.cos(C);
SI1 = Math.sin(A);
SI2 = Math.sin(B);
SI3 = Math.sin(C);
for (int I = J; I<N; I+=N1) {
I1 = I + N2;
I2 = I1 + N2;
I3 = I2 + N2;
R1 = X[I] + X[I2];
R3 = X[I] - X[I2];
S1 = Y[I] + Y[I2];
S3 = Y[I] - Y[I2];
R2 = X[I1] + X[I3];
R4 = X[I1] - X[I3];
S2 = Y[I1] + Y[I3];
S4 = Y[I1] - Y[I3];
X[I] = R1 + R2;
R2 = R1 - R2;
R1 = R3 - S4;
R3 = R3 + S4;
Y[I] = S1 + S2;
S2 = S1 - S2;
S1 = S3 + R4;
S3 = S3 - R4;
X[I1] = CO1*R3 + SI1*S3;
Y[I1] = CO1*S3 - SI1*R3;
X[I2] = CO2*R2 + SI2*S2;
Y[I2] = CO2*S2 - SI2*R2;
X[I3] = CO3*R1 + SI3*S1;
Y[I3] = CO3*S1 - SI3*R1;
}
}
}
// Radix-4 bit-reverse
double T;
int J = 0;
N2 = N>>2;
for (int I=0; I < N-1; I++) {
if (I < J) {
T = X[I];
X[I] = X[J];
X[J] = T;
T = Y[I];
Y[I] = Y[J];
Y[J] = T;
}
N1 = N2;
while ( J >= 3*N1 ) {
J -= 3*N1;
N1 >>= 2;
}
J += N1;
}
}
这是Sidney Burrus最初的Radix-4 DIF FORTRAN代码：
Radix-4，DIF，一只蝴蝶FFT
#2 楼
首先，您假定的“基数4蝴蝶”是4点DFT，而不是FFT。它具有16个复杂的（即：N平方）运算。典型的4点FFT仅具有Nlog（base 2）N（对于N = 4，= 8）。其次，您有一些假设的w [] .r和w [] .i“比例”因子不属于。也许您是从较大的图中显示的radix-4蝴蝶获得的。这样的蝴蝶将在其上附加一些级间旋转，但它们实际上并不是蝴蝶的一部分。当设计用于负指数FFT时，四点FFT仅具有-j的内部蝴蝶形。编译器；输出附加在代码末尾）：首先，打印输出输入数据（4个实数，4个虚数）。然后采用4点DFT。结果（yr []和yi []加上安培/相位）被打印出来。由于进行DFT时不会覆盖r []和i []原始数据，因此这些输入将重新用作4点FFT的输入。注意，后者比DFT具有更少的+/-操作。FFT的代码不是特别优雅也不高效-有许多种蝶形运算的方法。上面的代码对应于Rabiner和Gold的书“数字信号处理的理论和应用”（第580页，图10.9）中所示的四个基数为2的蝴蝶，并对其进行了修改以反映负指数（用于书中的数字是肯定的）。请注意，代码中的-j仅旋转一圈，并且不需要乘法（这是一个交换/符号更改）。它们与DFT相同
最后，将来自FFT的缩放结果用作逆FFT的输入。这可以通过“交换”或“反向列表”方法来完成（即：如果FFT（r，i）是正向FFT，则FFT（i，r）是逆向）–当然，只要能够处理复杂的输入/输出（换句话说），没有“仅真实的”例程（通常假定虚拟输入为零）。此方法已在25年前的
P中进行了描述。 Duhamel，B。Piron，J.M。Etcheto，“论逆DFT的计算”，《 IEEE Transactions on声学，语音和信号处理》，第1卷。 1988年2月36日，第285-286页。
然后将求反的结果打印出来。与原始输入数据相同。
评论
您还可以向我们展示每个示例的一些输入和输出数据示例吗？除了位反转顺序问题外，还有2倍或4倍的差异-一些实现对正向fft进行缩放，某些实现对反向fft进行缩放，而某些对这两种都进行缩放...
这不是一个重新排序的问题，因为据我了解，重新排序会置换y的条目。如果更改ang = -2 * PI，我可以解决此问题；而不是ang = -2 * PI /（fp）4;我不需要重新排列术语，我的一致性测试与dft传递的错误为0。我认为这相当于旋转因子90度相移。但是，这似乎与数学不一致……我想念的是什么？