疎信号に対する高速フーリエ変換
M2013SS006 溝口 竜二
指導教員：杉浦洋

1

はじめに

2.1

離散フーリエ変換 (DFT)

離散的時間の信号の離散的周波数スペクトルは離散フー高速フーリエ変換 (FFT) は基本的な数値アルゴリズム
リエ変換によって，計算量O(n log n) で求めることがでであり，n 次元信号の離散フーリエ変換 (DFT) は計算量2O(n ) から O(n log n) と高速で処理され，画像圧縮，信 きる．入力信号 x のフーリエスペクトルは以下のように号処理などの多くの応用分野で用いられている．信号が 定義される [1]．
特殊な条件を満たすときには，計算量をさらに減らすこ
1 ∑x̂i = √xj ω ij(1)とが出来る．今回はその条件として，信号のスパース性
nj∈Znに着目する．多くの応用分野で，信号のフーリエ係数のほとんどは 0 またはそれに等しいものである．このよう 離散フーリエ逆変換もまた，以下のように定義される．
な信号の性質をスパース性という．信号サイズが n，非零1 ∑周波数の K のとき，通常の計算では，これら K 個の周波xi = √xˆj ω −ij(2)n数成分を計算するのに，O(Kn) の計算量を必要とする．
j∈Zn
これより，計算量が小さいアルゴリズムを subliner time
√2πi/nは 1 の原始 n 乗根で，i = −1 である．
アルゴリズムといい，疎信号 FFT(SFFT) と呼ぶ．最初 ω = eの subliner time のアルゴリズムは，Hadamarld 変換に3 積和合同置換関してだが，[3] で提示された．すぐ後に，フーリエ変n を法として，奇数 σ ，整数 τ ∈ Zn が与えられたとき，換に対しての subliner time アルゴリズムが発見された
n次元ベクトルx に対する，積和合同置換 Pσ，τ は次のよ
[4][5][6][7][8]．
本研究では，MIT の Hassanieh らに提案された疎信号 うに定義される．
FFT(SFFT) のアルゴリズムを取り上げる [1]．この SFFT(Pσ，τ x)i = xσi+τ ．
(3)では，矩形関数とガウス関数を畳み込むことで得られるフラット窓関数をフィルタを用いる．このフラット窓関 置換のフーリエスペクトルは以下のように表せる．
数は，大きなフーリエ係数はパス領域内で区切られ，そ−τ i\(P．
(4)の周波数を特定する必要がある．そのサンプリングしたσ，
τ x)σi = x̂i ω
データに対して低次元 FFT(高速フーリエ変換) を行うことで，疎信号の重要でない部分とそうでない部分を容易に選り分けられるようにした．重要なデータを抽出すると，元信号の積和合同置換された位置データを抽出することができる．また，積和合同置換を行うことの利点は一つ信号データを指定した回数の異なる実行結果を得られることである．そこで得られた実行結果から，実行回数が半分以上出現した位置データの箇所だけを推定することで，元信号の計算量を減らし，元信号の近似解を推定する理論である. その領域だけで，すべての周波数を間接的に推定するアルゴリズムとして実行して，再現可能性や計算量の少ない再現アルゴリズムの有用性を確認す
る [2]．

2

変換関係と座標関係の表記

入力信号を x ∈ Cn とし，その記号と定義式を合わせて書く．Zn = {0，1，…，n−1} とし，x = (xj )n−1
j=0 は周期 n の∞無限列 (xj )j=−∞ の部分列と見なす．xj = xmod(j ，n) (j ∈Z) である．mod(j ，n) ∈ Zn は j を n で割った余りである．ベクトル x の台は supp(x) = {j ∈ Zn |xj ̸= 0} と定義する．

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.22004006008001000

図 1 変換前 (周波数領域)

4200400600800図 2 変換後 (周波数領域)

標準窓関数

本研究では，標準窓関数から派生するフラット窓関数を用いる [1]．デジタル信号処理では，通過域 ϵ，リップルの振幅の幅 δ ，窓の幅 w をパラメータとして持つ (ϵ，δ ，w)
標準窓関数 F ∈ Zn を以下のように定義する．
ω台 supp(F ) ⊆ [− ω2 ，2 ] をもち，対称 (F−j = Fj ) で，

• F̂0 = 1
• F̂i > 0i ∈ [−ϵn，ϵn]• |F̂i | < δi∈/ [−ϵn，ϵn]1000

フラット窓関数

57.1

Inner Loop

1 2，δ ，m) フ
任意の δ ，m = O(B log nδ ) に対して，( B
阻止域 ϵ，通過域 ϵ′ ，リップルの振幅の幅 δ ，窓の幅 w をB，′n
ラット窓関数を G とし，Location Loop は，パラメータパラメータとして持つ (ϵ，ϵ，δ ，w) フラット窓関数 G ∈ Z
ωd が与えられ，各非ゼロ係数を高い確率で含む dKnを以下のように定義する．台 supp(G) ⊆ [− ω2 ，B の集2 ] をもち，合 I ⊂ [n] を出力し，Estimation Loop は集合 I ⊂ [n] が対称 (G−j = Gj ) で，′′与えられ，x̂I は高い確率で推定する．
• Ĝi ∈ [1 − δ ，1 + δ]i ∈ [−ϵ n，ϵ n]• |Ĝi | < δ

1. 奇数 σ ∈ Zn と任意の整数 τ をランダムに選択する．

i∈/ [−ϵn，ϵn]2. 積和合同置換
入力ベクトル x を置換する x 7→ y (i ∈ Zn )−τ i\yi = (Pσ，τ x)i = xσi+τ ， ŷσi = (Pσ，τ x)σi = x̂i ω

2.0

3. フィルタの適用
フラット窓関数 G を用いて，置換されたベクトルを以下のように計算する．
zi = Gi · yi = Gi xσi+τ ，supp(z) ⊆ supp(G)

1.5

1.0

0.5

200

400600800

1000

4. サブサンプリング
置換されたフィルタをかけたベクトル z を B 個のビnンの集合に切り取る．L ≪ Bとして，
∑L
∑
n/B−1
∼
wl = l=−L zl+Bj =zl+Bj (l ∈ ZB )j=0B−1
5. 低次元√FFT：ŵl = (ŵ)l=0 を計算．
ŵl ∼= B ẑl( n ) となる

図 3 フラット窓関数 (時間領域)

n1.00.8

0.6

7. Location Loop(ρ)n−1i=0 の順列とし，ŵi の成分を絶対値の大きい順に並べたものを．

0.40.2|ŵρ(0) | ≥ |ŵρ(1) |… ≥ |ŵρ(n−1) |200400600800とする．このとき，J を以下のように定義する．
J = {ρ(0)，…，ρ(K ′ )} (K ′ = dK − 1)
J に属する番号の，ビンに属する座標 i の集合を，1000

図 4 フラット窓関数 (周波数領域)

6B

6. 関数とオフセット関数
ハッシュ関数を定義する．
hσ : Zn → ZB ，hσ (i) = round( σiBn )オフセット関数を定義する．
nnnoσ : Zn → [− 2B，2B ]，oσ (i) = σi − B hσ (i)I = {i ∈ Zn | hσ (i) ∈ J}

低次元 FFT

ベクトルを x ∈ Cn に対し，n の約数である自然数パラメータを B としたとき，l ∈ ZB で，∑とする．

8. Estimation Loopi ∈ I に対して，x̂i を以下のように推定する．

n/B−1ql =x̂′i =xl+Bjj=0で定義される，B 次元ベクトル q = {q0 , q1 , · · · , qB−1 } ∈CB のフーリエ変換を q̂ = {qˆ0 , qˆ1 , · · · , qB−1
ˆ } ∈ CB とする．このとき，√Bq̂l = x̂l(n/B)(l ∈ ZB ).
n7SFFT(Sparse Fast Fourier Transform)

SFFT は innner loop と outer loop の二つのアルゴリズムで構成され，以下のように記述する．

ŵhσ(i) ω τ iĜOσ(i)Outer Loop(個人選抜)√nK
d = O( 1ϵ ) をもち，パラメータ ϵ，δ ，B = O( ϵ logn )，7.2

δinner loop の実行回数を L = O(log n) とする．
1. Location loop の実行
L 回の Location loop を実行し，座標 I1 ，．
．
．
，IL の L個の集合を返す．
2. カウント
座標 i をを見つけるごとに出現数 si をカウントする．
si = |{r|i ∈ Ir }|3. 座標の保持
Location Loop で，si ≥ L/2 の半分発生した座標 iを保持する．
I ′ = {i ∈ I1 ∪ … ∪ IL | si ≥する．

ha = {16, 15, 10, 7}
WF = {1.00361, 0.914928, 0.914928, 0.999996}

L}．
2o = {5, 29, 29, 20}′′
4. フーリエ係数 x̂（i i ∈ I ）の推定′i ∈ I に対し，適当な inner loop で，Estimation Loop を 1 回実行し，フーリエ係数 x̂′i の推定を行う．

8SFFT の実行結果

8.1

Inner Loop の実行結果

estimatedx[[ii]] = {1. − 0.0850721I, 1., 1., 1.}

1.4
1.2

本研究では，n = 1024 の入力信号で，非零フーリエ係数 (K = 4) をもつ入力信号の離散フーリエ変換 x̂ の成分 (221, 445, 637, 956) = 1 で，その他は 0 というデータを作成し，B = 16 として，サブサンプリングと低次元 FFT を行い，SFFT を実行した．実行結果より，入力信号の離散フーリエ変換を高い確率で推定することができたが，アルゴリズムを実行していく中で，誤差が生じる結果の事例が生じた．Inner Loop は，location loop の回数を１回実行した結果をもとに，既知の位置データで，estimation loop を行い，入力信号を推定するアルゴリズムとする．

8.2

gF = {1., 0.914928, 0.914928, 0.999996}

1.0
0.8
0.6
0.4
0.2

200

非零フーリエ係数が，一つしか含まれないビンの Estimation Loop は成功する．次の例では，４つの非零フーリエ係数すべての推定が成功する．

WF = {1., 0.914928, 0.914928, 1.}

800

Outer Loop の結果

ランダム実験の成功例

L = 50 のとき

o = {5, 29, 29, 4}

位置データ = {176, 265, 377, 824}

gF = {1., 0.914928, 0.914928, 1.}

ha = {1, 11, 14, 6}

estimatedx[[ii]] = {1., 1., 1., 1.}

WF{1., 1., 1., 0.99696}
o = {17, 6, 12, 25}
gF = {1., 1., 1., 0.99696}

1.0

estimatedx[[iii]] = {1., 1., 1., 1.}
0.8

L = 100 のとき
0.6

位置データ = {64, 171, 199, 689}

ha = {12, 11, 10, 15}

0.4

WF{1., 1., 0.999996, 1}
0.2

o = {1, 8, 20, 17}
200

400

600

800

1000

gF = {1., 1., 0.999996, 1.}estimatedx[[iii]] = {1., 1., 1., 1.}

図 5 置換した信号とフラット窓関数を畳み込みした信号

8.3Estimation Loop の失敗例

非零フーリエ係数が複数含まれるビンは失敗する．次の例では，１番目の非零フーリエ係数の推定精度が劣化

1000

本研究では，location loop の実行回数 (L = 50, 100) で行い，実行回数の半分以上の位置データを抽出し，その位置データをもとに，estimation loop で推定した．

8.5

ha = {11, 2, 13, 7}600

図 6 置換した信号と窓関数を畳み込みした信号

8.4

Estimation Loop の成功例

400

ランダム実験の失敗例

入力信号の位置データ = {147, 338, 471, 537}

位置データが K 個正確に抽出される場合は，入力信号の離散フーリエ変換の非零成分のフーリエ係数を推定できるが，位置データが K 個以外の場合，正しく推定される部分もあるが，精度としてはあまり良くない結果となって推定された．

ha = {11, 12, 11, 6, 11}98.6

推定位置データ = {533, 146, 337, 470, 536}

まとめ（今後の課題）

WF = {0.408201, 0.980173, 0.408201, 0.999999, 0.408201}

SFFT のアルゴリズムの結果が上手くいかない事象を
o = {24, 27, 22, 19, 9}
もとに，誤って選ばれる周波数を抽出しないように，location loop の実行回数を増やすことが重要であると判明
gF = {0.998987, 0.980186, 0.999919, 0.999999, 1.}
した．極端に,location loop の実行回数を増やすと，計算estimated x[[iii]] = {0.124674 + 0.389131I, 0.999986 +
量が大きくなるので，実行結果を制限する必要がある．
0.000081384I, 0.0832468 − 0.399657I, 1., 0.0833951 +本研究のアルゴリズムは主に，location loop に依存し0.399592I}ている．入力信号に対して，擬似乱数をパラメータとした積和合同置換を実行しており，location loop 内の組織的8.7 誤って選ばれる周波数の期待値な理論が曖昧になっており，結果依存での location loop
L 回の Location loop のうち， L2 回以上の出現回数が を組むことで，理論的に追求する必要がある．フラットあったとき，非零周波数の位置データとして選択される． 窓関数と低次元 FFT は，計算量 (B log B) で出来る利点
この際に，固有値が 0 の周波数が誤って選ばれるときがあ もあるので，組織的に構築した位置データをもとに行う
る．その個数の平均値を計算すると，誤って選ばれる周波 事で、より効率のよい高速フーリエ変換のアルゴリズム数の期待値を，n：信号の大きさ，m：出現個数，B：ビン が出来るのではないかと考えられる．また，先行研究での個数，k：非ゼロフーリエ係数の個数，L：location loop ある MIT の研究者である Hassanieh らは，新たな改良しの出現回数，p：k がビンに含まれる確率とし，以下のよ た SFFT のアルゴリズムを提唱しているので，新しい理うな期待値で推定した．
論を追加するのも本研究のさらなる課題でもある．
()L∑L
E=(1 − p)(L−m) pm·n(5) 参考文献mLn= 2[1] H． Hassanieh,P． Indyk， D． Katabi,and E. Simple
and Practical Algorithm for Sparse Fourier Trans式 (5) を，L = 50, L = 100 のとき，半分の値 (25, 50) がform．SODA,2012
生じる期待値は，0.0823796,0.000172405 となり，L = 50
「圧縮センシングにおける完全再現十分条のときより L = 100 のとき，極端に半分の値の誤って選 [2] 和田山正，件について」, 日本神経回路学科誌，Vol． 17，No．
ばれる周波数の期待値が下がる．この結果，L = 100 回 の2(2010),63-69
実行結果でかなりの期待値で選択されないことがわかる．
[3] E. Kushilevitz and Y. Mansour. Learning decision8.8 location loop での見解
trees using the fourier spectrum. STOC, 1991
Inner Loop は，非零フーリエ係数の位置データを特定 [4] A. Gilbert, S. Guha, P. Indyk, M. Muthukrishnan,することができる．この結果より，location loop は非零and M. Strauss. Near-optimal sparse fourier repre周波数がかかる K 個のビンが重要であることがわかる．
sentations via sampling. STOC, 2002.
信号の離散フーリエ変換がフラット窓関数を畳み込む際[5] A. Akavia, S. Goldwasser, and S. Safra. Provに，非零成分の周波数を含むビンが重なる部分が生じ，位ing hard-core predicates using list decoding.FOCS,置データに誤差が生じる．このように，location loop は
pages146, 2003
１回では，失敗する可能性があるため，複数回する必要[6] A. Gilbert,M. Muthukrishnan,
and M.
がある．
Strauss.Improvedtimeboundsfornear-optimal
実行回数が L = 50 のとき，Outer loop の出現回数 sispace fourier representations.SPIE Conference,が実行回数の半分以上の位置データを推定する際に，誤っWavelets, 2005.て選ばれる周波数を推定してしまう場合が生じる．実行
回数が L = 100 とき，Outer loop を数十回実行しても， [7] M. A. Iwen. Combinatorial sublinear-time fourier al半分の 50 回に至ることはなかった.
gorithms.Foundations of Computational Mathematics, 10:303 338, 2010.
8.9 estimation loop での見解
[8] A. Akavia. Deterministic sparse fourier approxi非零成分の周波数を含むビンが K 個の場合，正確に非mation via fooling arithmetic progressions.COLT,零フーリエ係数を推定できる．また，ビンで衝突し，推pages381393, 2010.
定する箇所に大きな膨らみが生じた場合，位置データに誤差が生じることから，非零フーリエ係数は異なる解となり，正しく推定することが出来ないこともある．