DFT的計算

發文回覆瀏覽次數：2143

推到 Plurk!

推到 Facebook!

答題得分者是：GrandRURU

iamtonytu
一般會員

發表：1
回覆：0
積分：0
註冊：2015-06-16

發送簡訊給我

#1 引用回覆回覆

發表時間：2015-06-20 12:25:59

IP:180.217.xxx.xxx 訂閱

各位大大,小弟我爬過許多文章,想與各位確定DFT計算方式,在BCB裡執行是否有錯誤:
先不管計算量,小弟我怕自己邏輯錯誤,
1D的DFT

所以程式碼為

for (u=0;u {
bufferr=0.0;
bufferi=0.0
for (x=0;x {
angle=pi2*(u*x*xinvs);
bufferi=bufferi (array[x]*cos(angle));
bufferr=bufferr (array[x]*sin(angle));
}
F[u]=buffer;
}

在2D的DFT裡
double xinvs = 1.0 / X_Width ;
double yinvs = 1.0 / Y_Height;
for (u=0;u { for (v=0;v { buffer=0.0;
for (x=0;x {
for (y=0;y {
xangle=pi2*(u*x*xinvs);
yangle=pi2*(v*y*yinvs);
angle=xangle yangle;
buffer=buffer ((array[x][y])*(cos(angle)-sin(angle)));
}
}
F[u][v]=buffer;
this->Memo1->Lines->Add(AnsiString(F[u][v]));//之後顯示F[u][v] [Image1->Canvas->Pixels[x][y]=RGB(F[u][v],F[u][v],F[u][v]);]
}
}

找尋書說用2D計算方式,先算行(y axis)後再算列(x aixs).
但奇怪的是,
例如
double matrix[8]={1,0,-1,0,1,0,-1,0};
double M=8.0;
for (int u=1;u<8;u )
{ r=0.0;i=0.0;
for (int x=0;x<8;x )
{
r=r matrix[x]*cos(pi2*(double)x*(double)u/M);
i=i matrix[x]*sin(pi2*(double)x*(double)u/M);
}
this->Memo1->Lines->Add(AnsiString(r));
this->Memo2->Lines->Add(AnsiString(i));
}
書裡的答案是,0 , 0 , 0.5 , 0 , 0 , 0, 0.5 , 0
可是我算出來的答案卻是
r(real part):
7.26686506108409E-17
4
-3.29028254295038E-16
0
2.36174469737949E-15
4
9.34609377736817E-16
i (image part)
1.11022302462516E-16
-7.34763812293426E-16
4.44089209850063E-16
4.89842541528951E-16
6.66133814775094E-16
-2.20429143688028E-15
8.88178419700125E-16
到底是哪邊錯掉了,可否請各位大大幫忙,感謝

GrandRURU
站務副站長

發表：240
回覆：1680
積分：1874
註冊：2005-06-21

發送簡訊給我

#2 引用回覆回覆

發表時間：2023-04-28 08:47:05

IP:59.120.xxx.xxx 未訂閱

根據您提供的程式碼，似乎是 DFT 計算公式有些問題。以下是幾個可能的修正：

1D DFT 的計算公式應為：

[code cpp]
for (u=0;u {
F[u]=0.0;
for (x=0;x {
angle = -2.0 * pi * u * x / N;
F[u] = array[x] * complex(cos(angle), sin(angle));
}
}

[/code] 其中 complex 代表複數，因為 DFT 的結果是複數。

2D DFT 的計算公式應為：

[code cpp]
for (u=0;u {
for (v=0;v {
F[u][v] = complex(0.0, 0.0);
for (x=0;x {
for (y=0;y {
xangle = 2.0 * pi * u * x / N;
yangle = 2.0 * pi * v * y / N;
angle = xangle yangle;
F[u][v] = array[x][y] * complex(cos(angle), -sin(angle));
}
}
}
}

[/code] 其中 complex 代表複數，因為 DFT 的結果是複數。注意到這裡採用的是維納-庫因(Wiener-Khintchine)符號，也就是正負號相反。

由於浮點數計算的精度限制，最好不要直接比較浮點數是否相等。在顯示結果時，可以採用下面這種方式：

[code cpp]
if (fabs(F[u][v].real()) < 1e-6) F[u][v].real(0.0);
if (fabs(F[u][v].imag()) < 1e-6) F[u][v].imag(0.0);
this->Memo1->Lines->Add(AnsiString(F[u][v].real()));
this->Memo2->Lines->Add(AnsiString(F[u][v].imag()));

[/code] 其中 fabs 代表浮點數的絕對值，1e-6 是一個很小的數，相當於 $10^{-6}$。當一個複數的實部或虛部的絕對值小於這個數時，將其設為 0.0。
希望這些修正可以幫助您找到問題所在。

編輯記錄

GrandRURU 重新編輯於 2023-04-28 08:48:59，註解無‧
GrandRURU 重新編輯於 2023-04-28 08:50:50，註解無‧

系統時間：2024-04-25 9:22:07

聯絡我們 | Delphi K.Top討論版

本站聲明

1. 本論壇為無營利行為之開放平台，所有文章都是由網友自行張貼，如牽涉到法律糾紛一切與本站無關。
2. 假如網友發表之內容涉及侵權，而損及您的利益，請立即通知版主刪除。
3. 請勿批評中華民國元首及政府或批評各政黨，是藍是綠本站無權干涉，但這裡不是政治性論壇！