作者jim97260 (區欠足易)
看板C_and_CPP
標題[問題] DFT演算法的改寫
時間Tue May 21 17:08:14 2013
開發平台(Platform):(Ex: VC++, GCC, Linux, ...)
VC++
額外使用到的函數庫(Library Used):(Ex: OpenGL, ...)
No
問題(Question):
網路上找到的DFT演算法如下:
=============================================================================
for( x = 0; x < Size; ++x ) {
for ( int Step = 1; Step < Size; Step <<= 1 ) {
const int Jump = Step << 1;
const float delta = PI / float( Step );
const float Sine = sin( delta * .5 );
const complex Multiplier( -2. * Sine * Sine, sin( delta ));
complex Factor( 1. );
for ( int Group = 0; Group < Step; ++Group ) {
for ( int Pair = Group; Pair < Size; Pair += Jump ) {
const int Match = Pair + Step;
const complex Product( Factor * input[ Size * x + Match ] );
input[Size * x + Match] = input[Size * x + Pair] - Product;
input[Size * x + Pair] += Product;
}
Factor = Multiplier * Factor + Factor;
}
}
}
=============================================================================
大概說明:
Step是他每個處理的迴圈,是power of 2
Jump是每個Pair間的距離,也是power of 2
中間那些const是一些參數可以不用管
Group是迴圈裡的一個組,有幾個Group取決於Step,Factor會隨著Group而變化
ex:
------------------------
Size = 8
------------------------
Step 1, Jump = 2
Group 0 : Pair 0 (0, 1), Pair 1 (2, 3), Pair 2 (4, 5), Pair 3 (6, 7)
------------------------
Step 2, Jump = 4
Group 0 : Pair 0 (0, 2), Pair 1 (4, 6)
Group 1 : Pair 0 (1, 3), Pair 1 (5, 7)
------------------------
Step 4, Jump = 8
Group 0 : Pair 0 (0, 4)
Group 1 : Pair 0 (1, 5)
Group 2 : Pair 0 (2, 6)
Group 3 : Pair 0 (3, 7)
------------------------
在這種情況下,我如果要單獨計算某個位置是處於第幾個Pair或Group
以及他是Pair還是Match會很困難,因此我將這個演算法改寫成下面這樣:
===========================================================================
for( x = 0; x < Size; ++x ) {
for( y = 0; y < Size; ++y ) {
for ( Step = 1; Step < Size; Step <<= 1 ) {
int Jump = Step << 1;
// 判斷是Pair還是Match,Match的話不做動作
if (y % Jump < Step) {
float delta = PI / float( Step );
float Sine = sin( delta * .5 );
complex Multiplier( -2. * Sine * Sine, sin( delta ));
complex Factor( 1. );
// 判斷是第幾個Group
for ( int Group = 0; Group < ((y % Jump) - 1); ++Group ) {
Factor = Multiplier * Factor + Factor;
}
// Pair Match做運算
int Match = y + Step;
complex Product( Factor * input[ Size * x + Match ] );
input[ Size * x + Match ] = input[ Size * x + y ] - Product;
input[ Size * x + y ] += Product;
}
}
}
}
============================================================================
這樣整個演算法應該要變成下面這樣:
--------------------------------------------------
Size = 8
--------------------------------------------------
Step 1, Jump = 2
y = 0, y % Jump = 0 < Step, Group = 0, Match = 1
y = 1, y % Jump = 1 >= Step, Is Match, Do Nothing
y = 2, y % Jump = 0 < Step, Group = 0, Match = 3
y = 3, y % Jump = 1 >= Step, Is Match,Do Nothing
y = 4, y % Jump = 0 < Step, Group = 0, Match = 5
y = 5, y % Jump = 1 >= Step, Is Match,Do Nothing
y = 6, y % Jump = 0 < Step, Group = 0, Match = 7
y = 7, y % Jump = 1 >= Step, Is Match,Do Nothing
--------------------------------------------------
Step 2, Jump = 4
y = 0, y % Jump = 0 < Step, Group = 0, Match = 2
y = 1, y % Jump = 1 < Step, Group = 1, Match = 3
y = 2, y % Jump = 2 >= Step, Is Match,Do Nothing
y = 3, y % Jump = 3 >= Step, Is Match,Do Nothing
y = 4, y % Jump = 0 < Step, Group = 0, Match = 6
y = 5, y % Jump = 1 < Step, Group = 1, Match = 7
y = 6, y % Jump = 2 >= Step, Is Match,Do Nothing
y = 7, y % Jump = 3 >= Step, Is Match,Do Nothing
--------------------------------------------------
Step 4, Jump = 8
y = 0, y % Jump = 0 < Step, Group = 0, Match = 4
y = 1, y % Jump = 1 < Step, Group = 1, Match = 5
y = 2, y % Jump = 2 < Step, Group = 2, Match = 6
y = 3, y % Jump = 3 < Step, Group = 3, Match = 7
y = 4, y % Jump = 4 >= Step, Is Match,Do Nothing
y = 5, y % Jump = 5 >= Step, Is Match,Do Nothing
y = 6, y % Jump = 6 >= Step, Is Match,Do Nothing
y = 7, y % Jump = 7 >= Step, Is Match,Do Nothing
--------------------------------------------------
但是出來的結果卻不是正確的,所以想請問版上大大我的改寫是否有那裡出問題了嗎?
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 118.168.123.210
※ 編輯: jim97260 來自: 118.168.123.210 (05/21 17:14)
推 jackylu63:蝴蝶 05/21 22:12
推 hijamoya:CVDFT xddd 05/22 13:05
推 hexalift:心好像有點猛啊 05/27 18:40