重畳加算法のソースを表示

<!-- 英語版記事 [[:en:Overlap-add method]] 04:36, 10 December 2012 (UTC) (および同初版 と ドイツ語版記事 [[:de:Segmentierte Faltung]] 18:12, 2. Nov. 2012 (UTC)) 基づく 日本語版 -->
<!-- The '''overlap&ndash;add method (OA, OLA)''' is an efficient way to evaluate the discrete [[convolution]] of a very long signal <math>x[n]</math> with a [[finite impulse response]] (FIR) filter <math>h[n]</math>''':''' -->
'''重畳加算法''' (ちょうじょうかさんほう、{{lang-en|[[w:en:overlap-add method|OverLap-Add method]], OA, OLA}}) とは、非常に長い信号 <math>\mathbf{x}</math> と [[有限インパルス応答|FIR]]フィルタ <math>\mathbf{h}</math> の 離散[[畳み込み]] <math>\mathbf{h} * \mathbf{x}</math> を分割して(効率的に)処理する手法である。

<!-- :<math>
\begin{align}
y[n] = x[n] * h[n] \ \stackrel{\mathrm{def}}{=} \ \sum_{m=-\infty}^{\infty} h[m] \cdot x[n-m]
= \sum_{m=1}^{M} h[m] \cdot x[n-m],
\end{align}</math>
-->
:<math>
\begin{align}
y[n]	= (\mathbf{h} * \mathbf{x})[n]
	=  \sum_{m=0}^{M-1} h[m] \cdot x[n-m]
\end{align}</math>

<!-- where ''h''[''m''] = 0 for ''m'' outside the region [1, ''M'']. -->
ここで 信号 <math>x[n]</math>は <math>n=1,\ldots,N_x</math>、フィルタ <math>h[m]</math>は <math>m=0,\ldots,M-1</math> 以外で0、また<math>M<N_x</math>である。

<!-- The concept is to divide the problem into multiple convolutions of ''h''[''n''] with short segments of <math>x[n]</math>''':''' -->
重畳加算法の基本アイデアは、長い信号 <math>\mathbf{x}</math> を区間長<math>L</math>で区切り、複数の短い断片 <math>\mathbf{x}_k</math> と フィルタ <math>\mathbf{h}</math> に関する複数の<!-- ブロック -->畳み込みに[[分割統治法|分割]]して、{{trnote|適切なブロック長の}}[[高速フーリエ変換|FFT]]の積として効率的に計算することにある。

<!-- :<math>x_k[n]  \ \stackrel{\mathrm{def}}{=} 
\begin{cases}
x[n+kL] & n=1,2,\ldots,L\\
0 & \textrm{otherwise},
\end{cases}
</math>

where ''L'' is an arbitrary segment length.  Then''':'''

:<math>x[n] = \sum_{k} x_k[n-kL],\,</math>
-->
:<math>\begin{align}
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{align}</math>

<!-- and ''y''[''n''] can be written as a sum of short convolutions''':''' -->
離散畳み込み<math>(\mathbf{h} * \mathbf{x})[n]</math>は、各区間の畳み込みの総和で表せる:

<!-- :<math>
\begin{align}
y[n] = \left(\sum_{k} x_k[n-kL]\right) * h[n] &= \sum_{k} \left(x_k[n-kL]* h[n]\right)\\
&= \sum_{k} y_k[n-kL],
\end{align}
</math>
-->
:<math>\begin{align}
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{align}</math><!--
y_k[n]	& \stackrel{\mathrm{def}}{=}
	    (\mathbf{h} * \mathbf{x}_k)[n]
	  = \sum_{m=0}^{M-1} h[m] \cdot x_k[n-m] \\
	-->

<!-- where &nbsp;<math>y_k[n] \ \stackrel{\mathrm{def}}{=} \ x_k[n]*h[n]\,</math>&nbsp; is zero outside the region [1,&nbsp;''L''&nbsp;+&nbsp;''M''&nbsp;&minus;&nbsp;1]. &nbsp;And for any parameter &nbsp;<math>N\ge L+M-1,\,</math>&nbsp; it is equivalent to the <math>N\,</math>-point [[circular convolution]] of <math>x_k[n]\,</math> with <math>h[n]\,</math>&nbsp; in the region&nbsp;[1,&nbsp;''N'']. -->
ここで各区間の畳み込み <math>y_k[n] \ \stackrel{\mathrm{def}}{=} \ (\mathbf{h} * \mathbf{x}_k)[n] </math>は 領域 [1, ''L''+''M''-1] 以外で 0 であり、任意の<!-- パラメータ --> <math>N \ge (L+M-1)</math> に関し、領域 [1, N] における <math>\mathbf{x}_k</math> と <math>\mathbf{h}</math> の <math>N</math>点[[巡回畳み込み]]と等価である。

<!-- The advantage is that the circular convolution can be computed very efficiently as follows, according to the [[Discrete_Fourier_transform#Circular_convolution_theorem_and_cross-correlation_theorem|circular convolution theorem]]''':''' -->
重畳加算法のアドバンテージは、この巡回畳み込みが[[離散フーリエ変換#畳み込み定理と相互相関定理|畳み込み定理]]により、次のような FFTの積の逆FFT として効率的に計算できる事にある:

<!-- {{NumBlk|:|<math>y_k[n] = \textrm{IFFT}\left(\textrm{FFT}\left(x_k[n]\right)\cdot\textrm{FFT}\left(h[n]\right)\right)</math>|{{EquationRef|Eq.1}}}} -->
{{NumBlk|:|<math>y_k[n] = \textrm{IFFT}\left(\textrm{FFT}\left(\mathbf{x}_k\right)\cdot\textrm{FFT}\left(\mathbf{h}\right)\right)[n]</math>|{{EquationRef|式1}}}}

<!-- where FFT and IFFT refer to the fast Fourier transform and inverse
fast Fourier transform, respectively, evaluated over <math>N</math> discrete
points. -->
ここで <math>\textrm{FFT}</math>と<math>\textrm{IFFT}</math>は(ブロック長<math>N</math>の)[[高速フーリエ変換]]と逆高速フーリエ変換である。
<!-- ドイツ語版 [[:de:Segmentierte Faltung]] 18:12, 2. Nov. 2012 (UTC) より引用 -->FFTアルゴリズムによっては、(巡回畳み込み計算のために)FFTブロック長<math>N</math>を調整する事が理に適っている。例えば[[高速フーリエ変換#クーリー–テューキー型FFTアルゴリズム|Cooley-Tukey型FFTアルゴリズム]] (radix-2アルゴリズム) を使う場合、2の冪乗のブロック長を選ぶと有利である:

:<math>N = L + M - 1 = 2^p, \quad p \in \mathbb{N}</math>

== アルゴリズム ==

<!-- [[Image:Depiction of overlap-add algorithm.png|frame|none|Figure 1: the overlap&ndash;add method]] -->
[[画像:Depiction of overlap-add algorithm.png|frame|none|図1: 重畳加算法]]

<!-- Fig. 1 sketches the idea of the overlap&ndash;add method. The
signal <math>x[n]</math> is first partitioned into non-overlapping sequences,
then the [[discrete Fourier transform]]s of the sequences <math>y_k[n]</math>
are evaluated by multiplying the FFT of <math>x_k[n]</math> with the FFT of
<math>h[n]</math>. After recovering of <math>y_k[n]</math> by inverse FFT, the resulting
output signal is reconstructed by overlapping and adding the <math>y_k[n]</math>
as shown in the figure. The overlap arises from the fact that a linear
convolution is always longer than the original sequences. In the early days of development of the fast Fourier transform, <math>L</math> was often  chosen to be a power of 2 for efficiency, but further development has revealed efficient transforms for larger prime factorizations of L, reducing computational sensitivity to this parameter. A [[pseudocode]] of the algorithm is the
following: -->
図1に、重畳加算法のアイデアを示す。
# 信号<math>\mathbf{x}</math> (長さNx) を区間長<math>L</math>で区切り、オーバーラップの無いk個の断片の列<math>\mathbf{x}_k</math>を得る。
# 各区間について、断片 <math>\mathbf{x}_k</math> (長さ≦L) と フィルタ <math>\mathbf{h}</math> (長さM) のFFT結果の積をとり逆FFTして、区間毎の畳み込み結果 <math>\mathbf{y}_k</math> (長さ≦L+M-1) を得る。<br/>(なお畳み込み結果は、元信号(区間長L)と比べ (フィルタの長さ-1) だけ長くなり、オーバーラップが生じる)
# 最終的な出力信号は、図に示すように、各区間の結果<math>\mathbf{y}_k</math> のオーバーラップ(重畳)部を加算しながら連結して得られる。
区間長<math>L</math>は、FFT開発初期にはしばしば効率上の理由でFFTのブロック長<math>N</math>が2の冪乗をとるよう調整されたが、更なる研究開発の結果、Nを大きな素数で素因数分解する効率的変換方法が明らかにされ、このパラメータに対する計算上の敏感さは低減された。{{clarify|date=2013年1月}} ({{trnote|[[高速フーリエ変換|FFT]]のブロック長<math>N</math>が2の冪乗の場合の計算量は<math>O(N\log_2{N})</math>、一般にブロック長が <math>N=\Pi{n_i}</math> と素因数分解できる場合の計算量は<math>O(N\sum{n_i})</math>であり、ブロック長<math>N</math>をゼロ・パディングで調整して計算量を最適化できる。}})
アルゴリズムの[[疑似コード]]は以下の通り:

<!--    '''Algorithm 1''' (''OA for linear convolution'')
    Evaluate the best value of N and L
    H = FFT(h,N)       <span style="color:green;">(''zero-padded FFT'')</span>
    i = 1
    '''while''' i <= Nx
        il = min(i+L-1,Nx)
        yt = IFFT( FFT(x(i:il),N) * H, N)
        k  = min(i+N-1,Nx)
        y(i:k) = y(i:k) + yt    <span style="color:green;">(''add the overlapped output blocks'')</span>
        i = i+L
    '''end'''
-->
    '''アルゴリズム 1''' (''重畳加算法(OLA)による線形畳み込み'')
    使用するFFTアルゴリズムに応じてFFTブロック長 N と 分割区間長 L に最適値を設定
    H = FFT(h,N)       <span style="color:green;">(''ゼロ・パディングしたFFT'' )</span>
    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    <span style="color:green;">(''オーバーラップ区間を加算'' )</span>
        i = i+L
    '''end'''

== 巡回畳み込みへの応用 ==

<!-- When sequence ''x''[''n''] is periodic, and ''N''<sub>''x''</sub> is the period, then ''y''[''n''] is also periodic, with the same period. &nbsp;To compute one period of y[n], Algorithm 1 can first be used to convolve ''h''[''n''] with just one period of ''x''[''n'']. &nbsp;In the region ''M'' ≤ ''n'' ≤ ''N''<sub>''x''</sub>, &nbsp;the resultant ''y''[''n''] sequence is correct. &nbsp;And if the next ''M''&nbsp;&minus;&nbsp;1 values are added to the first ''M''&nbsp;&minus;&nbsp;1 values, then the region 1 ≤ ''n'' ≤ ''N''<sub>''x''</sub> will represent the desired convolution.  The modified pseudocode is''':''' -->
一般に信号<math>\mathbf{x}</math>が周期的でその周期が<math>N_x</math>の場合、畳み込み結果<math>y[n]</math>も周期的で同じ周期をとる。
# 1周期分の<math>y[n]</math>は、ちょうど1周期分の信号 <math>x[n]</math> (長さNx) と フィルタ<math>h[n]</math> (長さM) の畳み込みで得られる。ここでは前述のアルゴリズム 1を使う。<br/>畳み込み結果<math>y[n]</math>は、区間 [''M'', ''Nx''] で正しい。
# 先頭のM-1個(区間[1, ''M''-1])の値に、末尾のM-1個(区間[''Nx''+1, ''Nx''+''M''-1])の値を加算すれば、区間 [1, ''Nx''] が正しい畳み込み結果を表すようになる。
変更した疑似コードは以下の通り:{{vn|date=2013年1月}}

<!--    '''Algorithm 2''' (''OA for circular convolution'')
    Evaluate Algorithm 1
    y(1:M-1) = y(1:M-1) + y(Nx+1:Nx+M-1)
    y = y(1:Nx)
    '''end'''
-->
    '''アルゴリズム 2''' (''重畳加算法(OLA)による巡回畳み込み'')
    アルゴリズム 1 を評価
    y(1:M-1) = y(1:M-1) + y(Nx+1:Nx+M-1)
    y = y(1:Nx)
    '''end'''

<!--     Algorithm 2 (OA for circular convolution)
    \,
    Evaluate Algorithm 1
    ML = \left\lfloor \right.(NFFT-1)/L\left.\right\rfloor 
    i = Nx-L+1
    k = NFFT - L
    while k >= 1
        il = i+L-1
        yt = IFFT( FFT(x(i:il,NFFT)) * H )
        y(1:k) = y(1:k) + yt(NFFT-k+1:NFFT)
        k  = k-L
        i  = i-L
    end
-->
   【参考】[[:en:Special:PermaLink/137104117#Circular convolution with the Overlap-add method|英語版記事初版]]のアルゴリズム 2
    <span style="color:green;">(''アルゴリズム 1 の必要範囲を引用'' )</span>
        使用するFFTアルゴリズムに応じてFFTブロック長 N と 分割区間長 L に最適値を設定
        H = FFT(h,N)       <span style="color:green;">(''ゼロ・パディングしたFFT'' )</span>
    <span style="color:green;">(''ここまで引用'' )</span>
    ML = floor((N-1)/L)    <span style="color:green;">(''未使用'' )</span>
    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

== 計算量 ==
{{Notice|1=訳注: [[:en:Special:PermaLink/527298561#Cost of the overlap-add method|英語版記事 Overlap&ndash;add method 04:36, 10 December 2012 (UTC) 版のこの章]]にはいくつか問題があります:
<!-- * この章では、重畳加算法の目的は計算量削減よりも逐次分割処理のウェイトが高く、オーバーヘッドが充分小さい事を示せば充分に見える. -->
* 計算量<math>C_{OA}, C_s</math>が根拠不明 &mdash; 特に<math>C_s</math>の表式は、radix-2アルゴリズムの計算量に過ぎない.
* 計算量の比較意図が不明 &mdash; 重畳加算法のブロック畳み込みは巡回畳み込みと同じ形なのでFFTで高速化できるというストーリーなのに、重畳加算法を変形した巡回畳み込みAlgorithm 2と、高速化の鍵であるFFT処理(Eq.1)による(巡回畳み込み)アルゴリズム(詳細不明)を比較しようとしている.
* 評価基準の説明に不整合があり、"standard ..."の指す具体的アルゴリズムが不明.
** 本文記述: "''standard circular convolution using  Eq.1''" &mdash; Algorithm 2とは別の、Eq.1に基づくFFTベースの循環畳み込みの標準アルゴリズム
** Figure 2キャプション: "''time required by  Eq.1''" &mdash; Eq.1 (="Algorithm 1"相当): 循環畳み込みアルゴリズムではない
** [[:File:Gain oa method.png|図面]]の解説: "''standard FFT method''" &mdash; 標準FFT: 循環畳み込みアルゴリズムではない
** 本文中の計算量<math>C_s</math>: <math>N_x</math>が2の冪乗の場合の radix-2アルゴリズムの形 &mdash; 重畳加算法より計算量が小さいアルゴリズムが"standard"と呼ばれている
* 評価対象"Algorithm 2"の計測と図面作成後に、アルゴリズムが書き換えられている.
** [[:en:Special:PermaLink/137104117#Cost_of_the_Overlap-add_method|記事初版のFigure 2キャプション]]によれば、[[:en:Special:PermaLink/137104117#Circular_convolution_with_the_Overlap-add_method|20:11, 9 June 2007 版Algorithm 2]]の計算時間を実測(→図面に反映)
** [[:File:Gain oa method.png|図面]]は[[:c:File:Gain_oa_method.png#Metadata|02-Jun-2007 20:20:47]]作成
** その後、[[:en:Special:Diff/225200777#Circular_convolution_with_the_Overlap-add_method|12:42, 12 July 2008 (UTC)]]にAlgorithm 2が書き換えられている
}}

<!-- The cost of the convolution can be associated to the number of complex
multiplications involved in the operation. The major computational
effort is due to the FFT operation, which for a radix-2 algorithm
applied to a signal of length <math>N</math> roughly calls for <math>C=\frac{N}{2}\log_2 N</math>
complex multiplications. It turns out that the number of complex multiplications
of the overlap-add method are: -->
畳み込みの計算コストは、操作に関わる複素数乗算の回数と関連付ける事ができる。主要な計算量はFFT演算によるもので、Radix-2アルゴリズム({{trnote|ブロック長<math>N</math>が2の冪乗のアルゴリズム}})を長さ<math>N_x=N</math>の信号に適用する場合およそ <math>C = (N\log_2{N})/2</math>回の複素数乗算が行われる。重畳加算法における複素数乗算の回数は、FFTとフィルタ乗算とIFFTを考慮して:

:<math>C_{OA}=\left\lceil \frac{N_x}{N-M+1}\right\rceil
N\left(\log_2 N+1\right)\,</math>{{vn|date=2013年1月}}

<!-- <math>C_{OA}</math> accounts for the FFT+filter multiplication+IFFT operation. 
-->
<!-- The additional cost of the <math>M_L</math> sections involved in the circular
version of the overlap&ndash;add method is usually very small and can be
neglected for the sake of simplicity. The best value of <math>N</math>
can be found by numerical search of the minimum of <math>C_{OA}\left(N\right)=C_{OA}\left(2^m \right)</math>
by spanning the integer <math>m</math> in the range <math>\log_2\left(M\right)\le m\le\log_2 \left(N_x\right)</math>.
Being <math>N</math> a power of two, the FFTs of the overlap&ndash;add method
are computed efficiently. Once evaluated the value of <math>N</math> it
turns out that the optimal partitioning of <math>x[n]</math> has <math>L=N-M+1</math>.
For comparison, the cost of the standard circular convolution of <math>x[n]</math>
and <math>h[n]</math> is: -->
なお巡回行列版<!-- の重畳加算法 -->の<math>M_L</math>セクション ({{trnote|先頭と末尾の長さ(M-1)の区間？}}) の追加コストは、通常とても小さいので単純化のために無視できる。

FFTのブロック長<math>N</math>の最適値は、<math>\log_2{M}\le n\le \log_2{N_x}</math>の範囲で動かして<math>C_{OA}\left(N=2^n\right)</math>の最小値を整数<math>n</math>を(数値的に)探す事で得られる。<math>N</math>が2の冪乗であれば、<!-- 重畳加算法の -->FFTを効率的に計算できる。<math>N</math>の値が定まれば、信号<math>\mathbf{x}</math>を最適に区切る区間長<math>L=N-M+1</math>が定まる。
比較のため、普通の巡回畳み込みの計算量も示しておく:

:<math>C_S=N_x\left(\log_2 N_x+1\right)\,</math>{{vn|date=2013年1月}}

<!-- Hence the cost of the overlap&ndash;add method scales almost as <math>O\left(N_x\log_2 N\right)</math>
while the cost of the standard circular convolution method is almost
<math>O\left(N_x\log_2 N_x \right)</math>. However such functions accounts
only for the cost of the complex multiplications, regardless of the
other operations involved in the algorithm. A direct measure of the
computational time required by the algorithms is of much interest. 
Fig. 2 shows the ratio of the measured time to evaluate
a standard circular convolution using &nbsp;{{EquationNote|Eq.1}} with
the time elapsed by the same convolution using the overlap&ndash;add method
in the form of Alg 2, vs. the sequence and the filter length. Both algorithms have been implemented under [[Matlab]]. The
bold line represent the boundary of the region where the overlap&ndash;add
method is faster (ratio>1) than the standard circular convolution.
Note that the overlap&ndash;add method in the tested cases can be three
times faster than the standard method. -->
したがって重畳加算法の計算量はほぼ <math>O(N_x\log_2{N})</math>でスケールし、他方、普通の巡回畳み込みの計算量はほぼ <math>O(N_x\log_2{N_x})</math>である。({{trnote|<math>N>N_x</math>なので重畳加算法の方が計算量が多い事になってしまう。2つの計算量は要検証}}) しかしながら、この種の推計は複素数乗算の計算量だけ考慮しており、アルゴリズムに関わる他の処理は度外視している。各アルゴリズムに要する計算時間の直接測定こそ、より大きな関心事である。
図2に、{{EquationNote|式1}}による標準的な巡回畳み込み{{Clarify|date=2013年1月}}の計算時間と、''アルゴリズム 2''の形の重畳加算法による同様な畳み込みの計算時間の比を示す; 縦軸は信号長<math>N_x</math>の対数表示、横軸はフィルタ長<math>N_h=M</math>の対数表示で、計算時間の比を等高線で示している。両アルゴリズムは[[Matlab]]で実装した。青い太線は、重畳加算法の方が標準巡回畳み込みより速い(比>1)領域の境界線である。このケースでは重畳加算法は標準的手法より最大約3倍高速だった。{{vn|date=2013年1月}}

<!-- [[Image:gain oa method.png|frame|none|Figure 2: Ratio between the time required by &nbsp;{{EquationNote|Eq.1}} and the time required by the overlap&ndash;add Alg. 2 to evaluate
a complex circular convolution, vs the sequence length <math>N_x</math> and
the filter length <math>M</math>.]] -->
[[Image:gain oa method.png|frame|none|図2: {{EquationNote|式1}}の計算時間と、重畳加算法''アルゴリズム 2''の計算時間の比; 縦軸は信号長<math>N_x</math>の対数表示、横軸はフィルタ長<math>N_h=M</math>の対数表示]]

== 関連項目 ==

*{{仮リンク|重畳保留法}} ([[w:en:Overlap-save method|Overlap-save method]])

== 参考文献 ==

*{{Cite book
 | author    = Rabiner, Lawrence R.; Gold, Bernard
 | authorlink=
 | coauthors =
 | title     = Theory and application of digital signal processing
 | year      = 1975
 | publisher = Prentice-Hall
 | location  = Englewood Cliffs, N.J. 
 | isbn      = 0-13-914101-4
 | pages     = 63&ndash;67
}}
*{{Cite book
 | author    = Oppenheim, Alan V.; Schafer, Ronald W. 
 | authorlink= 
 | coauthors = 
 | title     = Digital signal processing 
 | year      = 1975 
 | publisher = Prentice-Hall 
 | location  = Englewood Cliffs, N.J.  
 | isbn      = 0-13-214635-5 
 | pages     =
}}
*{{Cite book
 | author    = Hayes, M. Horace 
 | authorlink= 
 | coauthors = 
 | title     = Digital Signal Processing
 | series    = Schaum's Outline Series
 | year      = 1999 
 | publisher = McGraw Hill 
 | location  = New York  
 | isbn      = 0-07-027389-8 
 | pages     =
}}

{{Reflist}}

== 外部リンク ==

{{DEFAULTSORT:ちようしようかさんほう}}
[[Category:信号処理]]
<!--[[Category:Transforms]]-->
[[Category:フーリエ解析]]
[[Category:数値解析]]
[[Category:数学に関する記事]]