高速フーリエ変換
をテンプレートにして作成
home
>
サイトマップ
開始行:
RIGHT:寄稿:東條遼平
* 制約と引替に高速処理 [#z53fd6ae]
画像をフーリエ変換する場合であれば、そこまで大きくない画...
まず通常の離散フーリエ変換の式を見てみます。
#ref(img13.png,nolink)
ここで
#ref(img34.png,nolink)
とし、一つ目の式を行列で表してみます。一般的に書けるので...
#ref(img35.png,nolink)
このとき
#ref(img36.png,nolink)
を使うと次のようになります。
#ref(img37.png,nolink)
ここで上の行列の真ん中で左右に分けて考えると、添字が偶数...
#ref(img38.png,nolink)
#ref(img39.png,nolink)
同じようすると最終的に次のようになります。
#ref(img40.png,nolink)
#ref(img41.png,nolink)
#ref(img42.png,nolink)
#ref(img43.png,nolink)
これを図で表すと次のようになります。
#ref(butterfly.png,nolink)
図1
図1を見てもわかる通りxの並びが変わっています。このように...
dfttimes = num;
power = power_two(num); //2の何乗かを返す
for(i=0; i<power-1; i++){
indexto = offset = 0;
while(indexto < num){
indexfrom = 0;
while(indexfrom < dfttimes){
rbuf[indexto] = re[indexfrom + offset];
ibuf[indexto] = im[indexfrom + offset];
indexto++;
indexfrom += 2;
if(indexfrom == dfttimes) indexfrom = 1;
}
offset += dfttimes;
}
for(j=0; j<num; j++){
re[j] = rbuf[j];
im[j] = ibuf[j];
}
dfttimes /= 2;
}
re, im, numは引数として渡されています。それぞれフーリエ変...
ここから少しソースがごちゃごちゃとするのですが、上の図1の...
unit_angle = 2.0 * PI / newnum;
dft = 2;
for(i=0; i<power; i++){
dftnum = newnum / dft;
step_angle = unit_angle * dftnum;
for(j=0; j<dftnum; j++){
angle = 0.0;
for(k=0; k<dft; k++){
indexto = j*dft + k;
if(k < dft/2){
indexfrom1 = indexto;
indexfrom2 = indexfrom1 + dft/2;
}else{
indexfrom2 = indexto;
indexfrom1 = indexto - dft/2;
}
rbuf = re[indexfrom2];
ibuf = im[indexfrom2];
temp_re[indexto] =
re[indexfrom1] + rbuf*cos(angle) + flag*ib...
temp_im[indexto] =
im[indexfrom1] + ibuf*cos(angle) - flag*rb...
angle += step_angle;
}
}
for(j=0; j<newnum; j++){
re[j] = temp_re[j];
im[j] = temp_im[j];
}
dft *= 2;
}
図1で矢印が交差している部分が視覚的に右から左へ3つのグル...
行列においてみられる対称性を利用し演算回数を減らしてきま...
- &ref(main.c);
- &ref(calculation.c);
- &ref(calculation.h);
終了行:
RIGHT:寄稿:東條遼平
* 制約と引替に高速処理 [#z53fd6ae]
画像をフーリエ変換する場合であれば、そこまで大きくない画...
まず通常の離散フーリエ変換の式を見てみます。
#ref(img13.png,nolink)
ここで
#ref(img34.png,nolink)
とし、一つ目の式を行列で表してみます。一般的に書けるので...
#ref(img35.png,nolink)
このとき
#ref(img36.png,nolink)
を使うと次のようになります。
#ref(img37.png,nolink)
ここで上の行列の真ん中で左右に分けて考えると、添字が偶数...
#ref(img38.png,nolink)
#ref(img39.png,nolink)
同じようすると最終的に次のようになります。
#ref(img40.png,nolink)
#ref(img41.png,nolink)
#ref(img42.png,nolink)
#ref(img43.png,nolink)
これを図で表すと次のようになります。
#ref(butterfly.png,nolink)
図1
図1を見てもわかる通りxの並びが変わっています。このように...
dfttimes = num;
power = power_two(num); //2の何乗かを返す
for(i=0; i<power-1; i++){
indexto = offset = 0;
while(indexto < num){
indexfrom = 0;
while(indexfrom < dfttimes){
rbuf[indexto] = re[indexfrom + offset];
ibuf[indexto] = im[indexfrom + offset];
indexto++;
indexfrom += 2;
if(indexfrom == dfttimes) indexfrom = 1;
}
offset += dfttimes;
}
for(j=0; j<num; j++){
re[j] = rbuf[j];
im[j] = ibuf[j];
}
dfttimes /= 2;
}
re, im, numは引数として渡されています。それぞれフーリエ変...
ここから少しソースがごちゃごちゃとするのですが、上の図1の...
unit_angle = 2.0 * PI / newnum;
dft = 2;
for(i=0; i<power; i++){
dftnum = newnum / dft;
step_angle = unit_angle * dftnum;
for(j=0; j<dftnum; j++){
angle = 0.0;
for(k=0; k<dft; k++){
indexto = j*dft + k;
if(k < dft/2){
indexfrom1 = indexto;
indexfrom2 = indexfrom1 + dft/2;
}else{
indexfrom2 = indexto;
indexfrom1 = indexto - dft/2;
}
rbuf = re[indexfrom2];
ibuf = im[indexfrom2];
temp_re[indexto] =
re[indexfrom1] + rbuf*cos(angle) + flag*ib...
temp_im[indexto] =
im[indexfrom1] + ibuf*cos(angle) - flag*rb...
angle += step_angle;
}
}
for(j=0; j<newnum; j++){
re[j] = temp_re[j];
im[j] = temp_im[j];
}
dft *= 2;
}
図1で矢印が交差している部分が視覚的に右から左へ3つのグル...
行列においてみられる対称性を利用し演算回数を減らしてきま...
- &ref(main.c);
- &ref(calculation.c);
- &ref(calculation.h);
ページ名:
home
>
Modified by
物理のかぎプロジェクト
PukiWiki 1.4.5_1
Copyright © 2001-2005
PukiWiki Developers Team
. License is
GPL
.
Based on "PukiWiki" 1.3 by
yu-ji
Powered by PHP 5.3.29HTML convert time to 0.002 sec.