重畳加算法 (ちょうじょうかさんほう、英語 : OverLap-Add method , OA, OLA
x
{\displaystyle \mathbf {x} }
FIR フィルタ
h
{\displaystyle \mathbf {h} }
畳み込み
h
∗
x
{\displaystyle \mathbf {h} *\mathbf {x} }
y
[
n
]
=
(
h
∗
x
)
[
n
]
=
∑
m
=
0
M
−
1
h
[
m
]
⋅
x
[
n
−
m
]
{\displaystyle {\begin{aligned}y[n]=(\mathbf {h} *\mathbf {x} )[n]=\sum _{m=0}^{M-1}h[m]\cdot x[n-m]\end{aligned}}}
ここで 信号
x
[
n
]
{\displaystyle x[n]}
n
=
1
,
…
,
N
x
{\displaystyle n=1,\ldots ,N_{x}}
h
[
m
]
{\displaystyle h[m]}
m
=
0
,
…
,
M
−
1
{\displaystyle m=0,\ldots ,M-1}
M
<
N
x
{\displaystyle M<N_{x}}
重畳加算法の基本アイデアは、長い信号
x
{\displaystyle \mathbf {x} }
L
{\displaystyle L}
x
k
{\displaystyle \mathbf {x} _{k}}
h
{\displaystyle \mathbf {h} }
分割 して、訳注 : 適切なブロック長のFFT の積として効率的に計算することにある。
x
k
[
n
]
=
d
e
f
{
x
[
n
+
k
L
]
,
n
=
1
,
2
,
…
,
L
0
,
otherwise
x
[
n
]
=
∑
k
x
k
[
n
−
k
L
]
,
{\displaystyle {\begin{aligned}x_{k}[n]\ &{\stackrel {\mathrm {def} }{=}}\ {\begin{cases}x[n+kL],&n=1,2,\ldots ,L\\0,&{\textrm {otherwise}}\end{cases}}\\x[n]&=\sum _{k}x_{k}[n-kL],\\\end{aligned}}}
離散畳み込み
(
h
∗
x
)
[
n
]
{\displaystyle (\mathbf {h} *\mathbf {x} )[n]}
y
[
n
]
=
(
h
∗
x
)
[
n
]
=
∑
m
=
0
M
−
1
(
h
[
m
]
⋅
∑
k
x
k
[
n
−
m
−
k
L
]
)
=
∑
k
(
h
∗
x
k
)
[
n
−
k
L
]
{\displaystyle {\begin{aligned}y[n]=(\mathbf {h} *\mathbf {x} )[n]&=\sum _{m=0}^{M-1}\left(h[m]\cdot \sum _{k}x_{k}[n-m-kL]\right)\\&=\sum _{k}(\mathbf {h} *\mathbf {x} _{k})[n-kL]\\\end{aligned}}}
ここで各区間の畳み込み
y
k
[
n
]
=
d
e
f
(
h
∗
x
k
)
[
n
]
{\displaystyle y_{k}[n]\ {\stackrel {\mathrm {def} }{=}}\ (\mathbf {h} *\mathbf {x} _{k})[n]}
L +M -1] 以外で 0 であり、任意の
N
≥
(
L
+
M
−
1
)
{\displaystyle N\geq (L+M-1)}
x
k
{\displaystyle \mathbf {x} _{k}}
h
{\displaystyle \mathbf {h} }
N
{\displaystyle N}
巡回畳み込み と等価である。
重畳加算法のアドバンテージは、この巡回畳み込みが畳み込み定理 により、次のような FFTの積の逆FFT として効率的に計算できる事にある:
y
k
[
n
]
=
IFFT
(
FFT
(
x
k
)
⋅
FFT
(
h
)
)
[
n
]
{\displaystyle y_{k}[n]={\textrm {IFFT}}\left({\textrm {FFT}}\left(\mathbf {x} _{k}\right)\cdot {\textrm {FFT}}\left(\mathbf {h} \right)\right)[n]}
(式1 )
ここで
FFT
{\displaystyle {\textrm {FFT}}}
IFFT
{\displaystyle {\textrm {IFFT}}}
N
{\displaystyle N}
高速フーリエ変換 と逆高速フーリエ変換である。
FFTアルゴリズムによっては、(巡回畳み込み計算のために)FFTブロック長
N
{\displaystyle N}
Cooley-Tukey型FFTアルゴリズム (radix-2アルゴリズム) を使う場合、2の冪乗のブロック長を選ぶと有利である:
N
=
L
+
M
−
1
=
2
p
,
p
∈
N
{\displaystyle N=L+M-1=2^{p},\quad p\in \mathbb {N} }
図1: 重畳加算法 図1に、重畳加算法のアイデアを示す。
信号
x
{\displaystyle \mathbf {x} }
L
{\displaystyle L}
x
k
{\displaystyle \mathbf {x} _{k}}
各区間について、断片
x
k
{\displaystyle \mathbf {x} _{k}}
h
{\displaystyle \mathbf {h} }
y
k
{\displaystyle \mathbf {y} _{k}}
最終的な出力信号は、図に示すように、各区間の結果
y
k
{\displaystyle \mathbf {y} _{k}}
区間長
L
{\displaystyle L}
N
{\displaystyle N}
[要説明 (訳注 : FFT のブロック長
N
{\displaystyle N}
O
(
N
log
2
N
)
{\displaystyle O(N\log _{2}{N})}
N
=
Π
n
i
{\displaystyle N=\Pi {n_{i}}}
O
(
N
∑
n
i
)
{\displaystyle O(N\sum {n_{i}})}
N
{\displaystyle N}
疑似コード は以下の通り:
アルゴリズム 1 (重畳加算法(OLA)による線形畳み込み )
使用するFFTアルゴリズムに応じてFFTブロック長 N と 分割区間長 L に最適値を設定
H = FFT(h,N) (ゼロ・パディングしたFFT )
i = 1
while i <= Nx
il = min(i+L-1,Nx)
yt = IFFT( FFT(x(i:il),N) * H, N)
k = min(i+M-1,Nx)
y(i:k) = y(i:k) + yt (オーバーラップ区間を加算 )
i = i+L
end
一般に信号
x
{\displaystyle \mathbf {x} }
N
x
{\displaystyle N_{x}}
y
[
n
]
{\displaystyle y[n]}
1周期分の
y
[
n
]
{\displaystyle y[n]}
x
[
n
]
{\displaystyle x[n]}
h
[
n
]
{\displaystyle h[n]}
y
[
n
]
{\displaystyle y[n]}
M , Nx ] で正しい。
先頭のM-1個(区間[1, M -1])の値に、末尾のM-1個(区間[Nx +1, Nx +M -1])の値を加算すれば、区間 [1, Nx ] が正しい畳み込み結果を表すようになる。 変更した疑似コードは以下の通り:[要検証 – ノート
アルゴリズム 2 (重畳加算法(OLA)による巡回畳み込み )
アルゴリズム 1 を評価
y(1:M-1) = y(1:M-1) + y(Nx+1:Nx+M-1)
y = y(1:Nx)
end
【参考】英語版記事初版 のアルゴリズム 2
(アルゴリズム 1 の必要範囲を引用 )
使用するFFTアルゴリズムに応じてFFTブロック長 N と 分割区間長 L に最適値を設定
H = FFT(h,N) (ゼロ・パディングしたFFT )
(ここまで引用 )
ML = floor((N-1)/L) (未使用 )
i = Nx-L+1
k = N - L
while k >= 1
il = i+L-1
yt = IFFT( FFT(x(i:il,N)) * H )
y(1:k) = y(1:k) + yt(N-k+1:N)
k = k-L
i = i-L
end
畳み込みの計算コストは、操作に関わる複素数乗算の回数と関連付ける事ができる。主要な計算量はFFT演算によるもので、Radix-2アルゴリズム(訳注 : ブロック長
N
{\displaystyle N}
N
x
=
N
{\displaystyle N_{x}=N}
C
=
(
N
log
2
N
)
/
2
{\displaystyle C=(N\log _{2}{N})/2}
C
O
A
=
⌈
N
x
N
−
M
+
1
⌉
N
(
log
2
N
+
1
)
{\displaystyle C_{OA}=\left\lceil {\frac {N_{x}}{N-M+1}}\right\rceil N\left(\log _{2}N+1\right)\,}
[要検証 – ノート なお巡回行列版の
M
L
{\displaystyle M_{L}}
訳注 : 先頭と末尾の長さ(M-1)の区間？
FFTのブロック長
N
{\displaystyle N}
log
2
M
≤
n
≤
log
2
N
x
{\displaystyle \log _{2}{M}\leq n\leq \log _{2}{N_{x}}}
C
O
A
(
N
=
2
n
)
{\displaystyle C_{OA}\left(N=2^{n}\right)}
n
{\displaystyle n}
N
{\displaystyle N}
N
{\displaystyle N}
x
{\displaystyle \mathbf {x} }
L
=
N
−
M
+
1
{\displaystyle L=N-M+1}
C
S
=
N
x
(
log
2
N
x
+
1
)
{\displaystyle C_{S}=N_{x}\left(\log _{2}N_{x}+1\right)\,}
[要検証 – ノート したがって重畳加算法の計算量はほぼ
O
(
N
x
log
2
N
)
{\displaystyle O(N_{x}\log _{2}{N})}
O
(
N
x
log
2
N
x
)
{\displaystyle O(N_{x}\log _{2}{N_{x}})}
訳注 :
N
>
N
x
{\displaystyle N>N_{x}}
式1 [要説明 の計算時間と、アルゴリズム 2 の形の重畳加算法による同様な畳み込みの計算時間の比を示す; 縦軸は信号長
N
x
{\displaystyle N_{x}}
N
h
=
M
{\displaystyle N_{h}=M}
Matlab で実装した。青い太線は、重畳加算法の方が標準巡回畳み込みより速い(比>1)領域の境界線である。このケースでは重畳加算法は標準的手法より最大約3倍高速だった。[要検証 – ノート
図2: 式1 アルゴリズム 2 の計算時間の比; 縦軸は信号長
N
x
{\displaystyle N_{x}}
N
h
=
M
{\displaystyle N_{h}=M}
Rabiner, Lawrence R.; Gold, Bernard (1975). Theory and application of digital signal processing . Englewood Cliffs, N.J.: Prentice-Hall. pp. 63–67. ISBN 0-13-914101-4 Oppenheim, Alan V.; Schafer, Ronald W. (1975). Digital signal processing . Englewood Cliffs, N.J.: Prentice-Hall. ISBN 0-13-214635-5 Hayes, M. Horace (1999). Digital Signal Processing . Schaum's Outline Series. New York: McGraw Hill. ISBN 0-07-027389-8
![{\displaystyle {\begin{aligned}y[n]=(\mathbf {h} *\mathbf {x} )[n]=\sum _{m=0}^{M-1}h[m]\cdot x[n-m]\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1e2c91b4c097802ebe36d43c05cac89f7e37b3de)
![{\displaystyle x[n]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/864cbbefbdcb55af4d9390911de1bf70167c4a3d)
![{\displaystyle h[m]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4bf03cff26a0e13bbd17cf88ee7b46b05c3df090)
![{\displaystyle {\begin{aligned}x_{k}[n]\ &{\stackrel {\mathrm {def} }{=}}\ {\begin{cases}x[n+kL],&n=1,2,\ldots ,L\\0,&{\textrm {otherwise}}\end{cases}}\\x[n]&=\sum _{k}x_{k}[n-kL],\\\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1f464f11275532584f82fb014da72829e2e87221)
![{\displaystyle (\mathbf {h} *\mathbf {x} )[n]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/5de9999be65eab67a1245ede733eaa6584609f7f)
![{\displaystyle {\begin{aligned}y[n]=(\mathbf {h} *\mathbf {x} )[n]&=\sum _{m=0}^{M-1}\left(h[m]\cdot \sum _{k}x_{k}[n-m-kL]\right)\\&=\sum _{k}(\mathbf {h} *\mathbf {x} _{k})[n-kL]\\\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/98ebea3ec4e9da47d8c04bd406c554c26deb6a9d)
![{\displaystyle y_{k}[n]\ {\stackrel {\mathrm {def} }{=}}\ (\mathbf {h} *\mathbf {x} _{k})[n]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8255e1bbf98795c87c954cae45c299f0734570e5)
![{\displaystyle y_{k}[n]={\textrm {IFFT}}\left({\textrm {FFT}}\left(\mathbf {x} _{k}\right)\cdot {\textrm {FFT}}\left(\mathbf {h} \right)\right)[n]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0dc07742e97e357e762f898bb560edf449253031)
![{\displaystyle y[n]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/305428e6d1fb59cd0163a7a96ace52292a262afa)
![{\displaystyle h[n]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/89981bbbb05ffd469eeadb828c18359965985e46)
