作者terrylau (國見比呂)
看板comm_and_RF
標題[討論] FFT C++ code
時間Thu Nov 16 22:19:47 2006
謝哪兩位得指教，我是有去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
1F：→ lovewa:我不曉得該怎麼回你，你應該先去看看FFT/IF 140.115.152.45 11/16 22:21
2F：→ lovewa:FT的數學表示法，然後再來看程式吧...=.= 140.115.152.45 11/16 22:22
3F：推 jammer:程式寫法也不難 218.184.208.98 11/16 22:25
4F：→ lovewa:因為上面的程式碼看起來只是在做那樣的事情 140.115.152.45 11/16 22:26
5F：→ lovewa:而且並不難，直覺他是直接把數學拆開來做.. 140.115.152.45 11/16 22:26
6F：推 kafai:沒有用recursive XD 140.112.245.86 11/30 12:38