[討論] FFT C++ code
謝哪兩位得指教，我是有去goodle過，就是看不懂，才來問！！
我手邊也有兩個版本的code,看不太懂
可以交一下嗎？？
typedef struct
{
double re;
double im;
}COMPLEX;
void fft(COMPLEX *x, int m)
{
COMPLEX w, s, t, u, v;
int n, l, le, le1, i, j, ip, nv2, nm1, k;
double rad;
n =(int) pow(2.0, m);
for (l = 1; l <= m; l++) {
le = (int)pow(2, m + 1 - l);
le1 = le / 2;
u.re = 1.0;
u.im = 0.0;
rad = PI / (float) le1;
w.re = cos(rad);
w.im = - sin(rad);
for (j = 0; j < le1; j++) {
for (i = j; i < n; i += le) {
ip = i + le1;
t.re = x[i].re + x[ip].re;
t.im = x[i].im + x[ip].im;
s.re = x[i].re - x[ip].re;
s.im = x[i].im - x[ip].im;
x[ip].re = s.re * u.re - s.im * u.im;
x[ip].im = s.re * u.im + s.im * u.re;
x[i] = t;
}
v.re = u.re * w.re - u.im * w.im;
v.im = u.re * w.im + u.im * w.re;
u = v;
}
}
/* bit reversal */
nv2 = n / 2;
nm1 = n - 1;
j = 0;
for (i = 0; i < nm1; i++) {
if (i < j) {
t = x[j];
x[j] = x[i];
x[i] = t;
}
k = nv2;
while (k < j + 1) {
j -= k;
k /= 2;
}
j += k;
}
for (i = 0; i < n; i++) {
x[i].re /= sqrt(n);
x[i].im /= sqrt(n);
}
return;
}
/*ifft*/
void ifft(COMPLEX *x, int m)
{
COMPLEX w, s, t, u, v;
int n, l, le, le1, i, j, ip, nv2, nm1, k;
double rad;
n =(int) pow(2.0, m);
for (l = 1; l <= m; l++) {
le = (int)pow(2, m + 1 - l);
le1 = le / 2;
u.re = 1.0;
u.im = 0.0;
rad = PI / (float) le1;
w.re = cos(rad);
w.im = sin(rad);
for (j = 0; j < le1; j++) {
for (i = j; i < n; i += le) {
ip = i + le1;
t.re = x[i].re + x[ip].re;
t.im = x[i].im + x[ip].im;
s.re = x[i].re - x[ip].re;
s.im = x[i].im - x[ip].im;
x[ip].re = s.re * u.re - s.im * u.im;
x[ip].im = s.re * u.im + s.im * u.re;
x[i] = t;
}
v.re = u.re * w.re - u.im * w.im;
v.im = u.re * w.im + u.im * w.re;
u = v;
}
}
/* bit reversal */
nv2 = n / 2;
nm1 = n - 1;
j = 0;
for (i = 0; i < nm1; i++) {
if (i < j) {
t = x[j];
x[j] = x[i];
x[i] = t;
}
k = nv2;
while (k < j + 1) {
j -= k;
k /= 2;
}
j += k;
}
for (i = 0; i < n; i++) {
x[i].re /= sqrt(n);
x[i].im /= sqrt(n);
}
return;
}
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 163.22.18.89
→
11/16 22:21, , 1F
11/16 22:21, 1F
→
11/16 22:22, , 2F
11/16 22:22, 2F
推
11/16 22:25, , 3F
11/16 22:25, 3F
→
11/16 22:26, , 4F
11/16 22:26, 4F
→
11/16 22:26, , 5F
11/16 22:26, 5F
推
11/30 12:38, , 6F
11/30 12:38, 6F