超急！ 誰能幫我解釋C程式，有重賞 ” 20 點 ”！

Question

★解釋越詳細越好！！

FFT是快速傅立葉轉換
Complex是複數(實部┼虛部)

主程式：

main()
{
COMPLEX y[8]={1000,0,1000,0,1000,0,1000,0,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}; /*rectangular input*/
int i, n = 8;
FFT(y,n); /*calls generic FFT function*/
for (i = 0; i {
*out_addr++ = (y[i]).real; /*real output component */
*out_addr++ = (y[i]).imag; /*imaginary output component*/
}
}

副程式：

/*FFT.C - FFT RADIX-2 USING DIF. FOR UP TO 512 POINTS */
#include "complex.h"
#include "twiddle.h"
void FFT(COMPLEX *Y, int N)
{
COMPLEX temp1,temp2;
int i,j,k;
int upper_leg, lower_leg;
int leg_diff;
int num_stages=0;
int index, step;
i=1;
do
{
num_stages+=1;
i=i*2;
} while (i!=N);

leg_diff=N/2;
step=512/N;
for (i=0;i {
index=0;
for (j=0;j {
for (upper_leg=j;upper_leg {
lower_leg=upper_leg+leg_diff;
temp1.real=(Y[upper_leg]).real + (Y[lower_leg]).real;
temp1.imag=(Y[upper_leg]).imag + (Y[lower_leg]).imag;
temp2.real=(Y[upper_leg]).real - (Y[lower_leg]).real;
temp2.imag=(Y[upper_leg]).imag - (Y[lower_leg]).imag;
(Y[lower_leg]).real=temp2.real*(w[index]).real-temp2.imag*(w[index]).imag;
(Y[lower_leg]).imag=temp2.real*(w[index]).imag+temp2.imag*(w[index]).real;
(Y[upper_leg]).real=temp1.real;
(Y[upper_leg]).imag=temp1.imag;
}
index+=step;
}
leg_diff=leg_diff/2;
step*=2;
}
j=0;
for (i=1;i<(N-1);i++) /*bit reversal for resequencing data*/
{
k=N/2;
while (k<=j)
{
j=j-k;
k=k/2;
}
j=j+k;
if (i {
temp1.real=(Y[j]).real;
temp1.imag=(Y[j]).imag;
(Y[j]).real=(Y[i]).real;
(Y[j]).imag=(Y[i]).imag;
(Y[i]).real=temp1.real;
(Y[i]).imag=temp1.imag;
}
}
return;
}

這是主程式main()前的部分，不小心漏掉了！

/*FFT8C.C - 8-POINT COMPLEX FFT PROGRAM. CALLS FFT.C */
#include "complex.h" /*complex structure definition */
extern void FFT(); /*FFT function */
volatile int *out_addr=(volatile int *)0x809802; /*out addr*/

副程式後段的蜂巢式迴圈是"位元反轉"(bit reversal)
但是我看不懂它的運作....

Answer 1

我回信給你 12小時 10 分了，你都沒回應！
看來你是不急了！

我要去忙我的了。
bye!

祝你好運。

2007-07-09 10:33:47 補充：
1. 你學長的碼不完整：少了 complex.h 和 twiddle.h
2. 你學長的碼有潛在性 bug：FFT 裡的第一個 do while 可能會 loop forever
3. 你學長的碼有多／錯的部分：那個 loop 看來是多餘的
4. 你學長對 C 應該很不熟！

2007-07-09 10:36:30 補充：
5. complex.h 我猜是 structCOMPLEX{doublereal;doubleimag;};
　我還可以補上去，twiddle.h 我就沒辦法了！
　它應該是在定義／計算 w[index]
　我 trace 了一下，index 可達 256 還是 512！
　w 是 FFT 的重點！
　更是硬體裡 FFT 加速的精神所在！
　沒有 w 的上限值，很難幫你弄出來！

2007-07-09 10:40:27 補充：
　（w 是 root of unity，不難算。
　　但要知道它的量，程式才會好寫、快！
　　這些資訊都沒有，總不致於要別人免費都幫你把這些找出來吧！）
6. FFT 的確可以用在整數上！
　但你沒說你為何要用 FFT，程式又不能跑！
　不知答案是否為 int！
　所以，最後的 volatile int *out_addr=(volatile int *)0x809802; 〝可能〞有錯！
　至少， volatile 應不應該這樣用，都是個大問題！
　就算能，809802 那位址可能是特殊硬體上的吧！？
　你這是要做什麼硬體也沒說！

2007-07-09 10:51:35 補充：
FFT 應用太廣了！
第一時間回信問你，你去上班等了超過 12 小時才回信。
你又沒說你要用在哪，程式又不能 run。
抱歉，幫不了你了！

這種不能完成的東西，到現在我都是經常熬整夜去征服它！
年輕時 (20年前)，我還可以 72 小時以上
　不睡、
　除了上廁所以外，不離開椅子（吃喝是別人端來／拿走）
去把它完成。
不是要你學我熬夜，那太傷身。
也不是要你像我以前一樣不懂事：虐待你身邊最親的人！
而是要你知道別人是怎麼在拚，才〝有希望〞真的學會。

下次要急著找人幫忙，不12小時以上沒消息！
萬一像這次一樣缺一堆，大概很難找到人幫你！

做過 數次 FFT、想幫你實在沒辦法的人 留

2007-07-09 10:54:55 補充：
你若要自己 K 懂，C.L.R.S. 的<>裡有一個相當不錯的 FFT in place recursion 的 algorithm！
因為它是 in place，很快！
〝看起來〞免bit reversal（神奇吧！）

你自己翻翻

2007-07-09 15:24:48 補充：
都沒人回答，你又來信要我不必客氣。
那我臉皮就厚一點。

回答內容如意見所述。

另外，若它可以跑，還方便 check 它的 bit reverse 對不對，
也比較有動力去了解。
萬一那 code 是錯的，看了半天．．．

所以，你自己加油囉！

有關FFT的程式，個人看過、認為講解最清淅、程式好懂、易寫出來的，就是意見裡寫的那本 C. L. R. S. 的書了。
去圖書館找找，拿來看看，或許，你自己可以從頭寫個 FFT，不用學長的。

加油！

Answer 2

完蛋了，在學校學的傅立葉公式完全忘光光…。

＝＿＝|||