反転M系列によるインパルス応答測定
Impulse Response Measurement Using Time-Reversed MLS
本稿は2018年8月25日にNHK放送技術研究所で開催された音楽音響研究会（日本音響学会の分科会のひとつ）で発表した内容に基づく。 証明に DFT を使用しているが、実際の装置設計や測定に DFT は必要なく、通常の FIR 演算のみでインパルス応答を測定できる。
はじめに
線形時不変系に $\pm 1$ に2値化したM系列を入力し、出力との相互相関を Hadamard 変換を用いて効率的に求め、スケーリング補正してインパルス応答を得る方法が知られている (Borish ら, 金田)。
ここでは最初にM系列のゲインとオフセットを補正して、DFT による振幅スペクトルが全周波数で1になるようにした周期信号を入力に用い、その時間軸を反転させた信号1周期との線形畳み込み（巡回的な相互相関と等価）によってインパルス応答を測定する方法について述べる。
定義
単位パルスを次のように定義する。
$$ \delta(n) = \left\{ \begin{array}{cl} 1, & n=0 \\ 0, & \text{その他} \end{array} \right. \tag{1} $$$N$ 点 DFT と逆 DFT の定義として以下を採用する。
$$ \begin{aligned} \text{DFT} &: X(k) = \sum_{n=0}^{N-1} x(n)\exp\!\left(-2\pi i \frac{nk}{N}\right) \\[6pt] \text{逆 DFT} &: x(n) = \frac{1}{N}\sum_{k=0}^{N-1} X(k)\exp\!\left(2\pi i \frac{nk}{N}\right) \end{aligned} \tag{2} $$周期 $N$ の周期関数の自己相関を $0\leq n\lt N$ の範囲で次のように定義する。
$$ r(n) = \sum_{i=0}^{N-1} x(i)\, x(i+n) \tag{3} $$全周波数で振幅スペクトルが1の周期信号を作る
ゲイン調整
周期 $N$ の2進M系列 $m(n)\in\{0,1\}$ から、次のように2値 $\pm\displaystyle\frac{1}{\sqrt{N}}$ を取る周期波形
$$ x(n) = \frac{1}{\sqrt{N}}\left\{2\, m(n)-1\right\} \tag{4}\label{xn} $$を作ると、その自己相関は次のように表せる (柏木)。
$$ r(n) = \left\{ \begin{array}{cl} 1, & n=0 \\ -\displaystyle\frac{1}{N}, & \text{その他} \end{array} \right. = \left(1+\frac{1}{N}\right)\delta(n)-\frac{1}{N} \tag{5}\label{rn} $$自己相関の DFT により $x(n)$ のパワースペクトルは、自己相関が偶関数であることから
$$ \begin{aligned} |\hat{X}(k)|^2 &= \sum_{n=0}^{N-1} r(n) \exp\!\left(-2\pi i \frac{nk}{N}\right) \\ &= \sum_{n=0}^{N-1} r(n) \cos\!\left(-2\pi \frac{nk}{N}\right) \\ &= \sum_{n=0}^{N-1} \left\{ \left(1+\frac{1}{N}\right)\delta(n)-\frac{1}{N} \right\} \cos\!\left(-2\pi \frac{nk}{N}\right) \quad (\text{式}\eqref{rn}\text{より}) \\ &= \left(1+\frac{1}{N}\right)\sum_{n=0}^{N-1} \delta(n) \cos\!\left(-2\pi \frac{nk}{N}\right) - \frac{1}{N} \sum_{n=0}^{N-1} \cos\!\left(-2\pi \frac{nk}{N}\right) \\ &= \frac{N+1}{N} - \delta(k) \end{aligned} \tag{6} $$$x(n)$ の振幅スペクトルはパワースペクトルの平方根であり、直流を除き一定値となる。
$$ |\hat{X}(k)| = \left\{ \begin{array}{ll} \displaystyle\frac{1}{\sqrt{N}}, & k=0 \\[8pt] \displaystyle\sqrt{\frac{N+1}{N}}, & 0\lt k\leq \displaystyle\frac{N-1}{2} \end{array} \right. = \left(\frac{1}{\sqrt{N}}-\sqrt{\frac{N+1}{N}}\right)\delta(k) + \sqrt{\frac{N+1}{N}} \tag{7}\label{AbsX} $$$x(n)$ に $\displaystyle\sqrt{\frac{N}{N+1}}$ を掛けた波形の振幅スペクトルは、式$\eqref{AbsX}$より
$$ \sqrt{\frac{N}{N+1}}\,|\hat{X}(k)| = \left(\frac{1}{\sqrt{N+1}}-1\right)\delta(k) + 1 \tag{8}\label{RootNN1X} $$オフセット調整
式$\eqref{RootNN1X}$の右辺第1項を移項したスペクトルを $Y(k)$ とすると、その振幅は全周波数で1になる。
$$ |Y(k)| = \sqrt{\frac{N}{N+1}}\,|\hat{X}(k)| + \left(1-\frac{1}{\sqrt{N+1}}\right)\delta(k) = 1 \tag{9}\label{all1} $$ここで $\left(1-\displaystyle\frac{1}{\sqrt{N+1}}\right)\delta(k)$ の逆 DFT は
$$ \frac{1}{N} \left(1-\frac{1}{\sqrt{N+1}}\right) \sum_{k=0}^{N-1} \delta(k) \exp\!\left(2\pi i \frac{nk}{N}\right) = \frac{1}{N} \left(1-\frac{1}{\sqrt{N+1}}\right) \tag{10} $$だから、$Y$ の逆 DFT は次のようになる。
$$ \begin{aligned} y(n) &= \sqrt{\frac{N}{N+1}}\,x(n) + \frac{1}{N}\left(1-\frac{1}{\sqrt{N+1}}\right) \\ &= \frac{\sqrt{N}}{\sqrt{N+1}}\cdot\frac{1}{\sqrt{N}}\left\{2\, m(n)-1\right\} + \frac{1}{N}\left(1-\frac{1}{\sqrt{N+1}}\right) \quad (\text{式}\eqref{xn}\text{より}) \\ &= \frac{2}{\sqrt{N+1}}\,m(n)+ \frac{1-\sqrt{N+1}}{N} \end{aligned} \tag{11} $$$y(n)$ の時間反転 $y(-n)$ の DFT は
$$ \sum_{n=0}^{N-1} y(-n)\exp\!\left(-2\pi i \frac{nk}{N}\right) = Y^*(k) \tag{12} $$であり、$y(n),\ y(-n)$ 各1周期に対する巡回畳み込みを $*$ で表すと以下が成り立つ。
$$ \begin{aligned} y(n) * y(-n) &= \frac{1}{N}\sum_{k=0}^{N-1} \hat{Y}(k)\hat{Y}^*(k) \exp\!\left(2\pi i \frac{nk}{N}\right) \\ &= \frac{1}{N}\sum_{k=0}^{N-1} |\hat{Y}(k)|^2 \exp\!\left(2\pi i \frac{nk}{N}\right) \\ &= \frac{1}{N}\sum_{k=0}^{N-1} \exp\!\left(2\pi i \frac{nk}{N}\right) \quad (\text{式}\eqref{all1}\text{より } |\hat{Y}(k)|^2=1) \\ &= \delta(n) \end{aligned} \tag{13}\label{cynyn} $$反転M系列を用いたインパルス応答測定
測定に使用する信号
$y(n),\ y(-n)$ から1周期を切り出した以下の孤立波形を考える。
$$ y_+(n) = \left\{ \begin{array}{ll} y(+n), & 0\leq n\lt N \\ 0, & \text{その他} \end{array} \right. \tag{14} $$ $$ y_-(n) = \left\{ \begin{array}{ll} y(-n), & 0\leq n\lt N \\ 0, & \text{その他} \end{array} \right. \tag{15} $$$y(n)$ は $y_+(n)$ を使って
$$ y(n) = \sum_{m=-\infty}^{\infty} y_+(n) * \delta(n-mN) \tag{16} $$と書けるので、これに $y_-(n)$ を畳み込むと
$$ \begin{aligned} y(n)*y_-(n) &= \left\{\sum_{m=-\infty}^{\infty} y_+(n)*\delta(n-mN)\right\}*y_-(n) \\ &= y_+(n)*y_-(n) * \sum_{m=-\infty}^{\infty} \delta(n-mN) \end{aligned} \tag{17} $$となり、孤立波形 $y_+(n)*y_-(n)$ と周期波形 $\displaystyle\sum_{m=-\infty}^{\infty} \delta(n-mN)$ の線形畳み込みは、周期 $N$ の巡回畳み込みになるので、式$\eqref{cynyn}$を適用できて
$$ y(n)*y_-(n) = \delta(n) * \sum_{m=-\infty}^{\infty} \delta(n-mN) = \sum_{m=-\infty}^{\infty} \delta(n-mN) \tag{18}\label{ynyn} $$と周期 $N$ の単位パルス列になる。
インパルス応答の測定原理
未知の線形時不変系のインパルス応答 $h(n)$ を測定することを考える。ただし、インパルス応答の長さは $N$ サンプル以下と仮定する。
$$ h(n) = 0,\quad n\lt 0 \text{ または } N\leq n \tag{19}\label{hncond} $$この系に周期信号 $y(n)$ を入力し、ノイズ $w(n)$ が重畳した出力 $f(n)$ が得られたとする。
$$ f(n) = h(n)*y(n) + w(n) \tag{20} $$$f(n)$ にゲインとオフセットを調整した反転M系列の1周期分 $y_-(n)$ を畳み込んだものを $g(n)$ とする。
$$ \begin{aligned} g(n) &= f(n)*y_-(n) \\ &= \left\{h(n)*y(n)+w(n)\right\}*y_-(n) \\ &= h(n)*y(n)*y_-(n) + w(n)*y_-(n) \\ &= \left\{h(n)*\sum_{m=-\infty}^{\infty} \delta(n-mN)\right\} + \underbrace{w(n)*y_-(n)}_{=w'(n)\text{と置く}} \quad (\text{式}\eqref{ynyn}\text{より}) \\ &= \left\{\sum_{m=-\infty}^{\infty} h(n-mN)\right\} + w'(n) \end{aligned} \tag{21}\label{fnyn} $$$g(n)$ を $N$ サンプル毎に区切って $J$ 回アベレージングしたものを $\hat{h}(n),\ 0\leq n\lt N$ とすると
$$ \begin{aligned} \hat{h}(n) &= \frac{1}{J}\sum_{j=0}^{J-1} g(n+jN) \\ &= \frac{1}{J}\sum_{j=0}^{J-1} \left[ \left\{\sum_{m=-\infty}^{\infty} h(n+jN-mN)\right\} + w'(n+jN) \right] \quad (\text{式}\eqref{fnyn}\text{より}) \end{aligned} \tag{22} $$条件式$\eqref{hncond}$が成り立ち、$0\leq n\lt N$ だから、$\displaystyle\sum_{m=-\infty}^{\infty}$ は $m=j$ の項だけが残り
$$ \hat{h}(n) = \frac{1}{J}\sum_{j=0}^{J-1} h(n) + \frac{1}{J}\sum_{j=0}^{J-1}w'(n+jN) \tag{23} $$系を線形時不変と仮定しているから、第1項の $h(n)$ は $j$ によらず確定で、$E[w'(n)]=0$ なら $J\to\infty$ で第2項は0になるので、以下が成り立つ。
$$ \lim_{J\to\infty}\hat{h}(n) = h(n) \tag{24} $$数値実験
計算法の正しさを確認するため、PC上の数値実験で、ノイズを加えずに次のインパルス応答 $h(n)$ を持つ系を推定する。
$$ h(n)=\{1.0,\ 0.9,\ 0.8,\ 0.7,\ 0.6,\ 0.5,\ 0.4,\ 0.3,\ 0.2,\ 0.1\} $$実験条件
インパルス応答の長さ10よりも長い $N=15$ のM系列を使い、$N$ が小さいので畳み込みには FIR を使い、計算はすべて倍精度浮動小数点で行う。
実験結果
図2にシミュレーション構成を、図3にシミュレーション結果を示す。
$g(n)$ の先頭2周期は乱れているが、3周期目以降には正しいインパルス応答 $h(n)$ が繰り返し得られている（誤差は最大で $2.22\times 10^{-16}$ 程度）。
考察
$g(n)$ の先頭2周期が乱れるのは、測定開始 $n=0$ 以前の系への入力が0であり、$n\lt 0$ に対しても入力の周期性を前提にしている計算原理と整合しないためである。
$N$ サンプルの信号同士を畳み込むと $2N-1$ サンプルになるので、先頭2周期では巡回畳み込みが成立せず、乱れを生じる。
このためアベレージングを行う場合は、先頭2周期（$2N$ サンプル）を除外する必要がある。
まとめ
M系列 $m(n)\in\{0,1\}$ を元にして、振幅スペクトルが全周波数で1になるようにゲインとオフセットを調整した2値の周期信号
$$ y(n) = \frac{2}{\sqrt{N+1}}\,m(n)+ \frac{1-\sqrt{N+1}}{N} $$を測定したい系に入力し、$y(n)$ の時間軸を反転させた1周期
$$ y_-(n) = \left\{ \begin{array}{ll} y(-n), & 0\leq n\lt N \\ 0, & \text{その他} \end{array} \right. $$を系の出力信号に畳み込むことで、系のインパルス応答が周期 $N$ で繰り返し得られることを導き、数値実験により計算の正しさを確認した。
付記
今回の実験はインパルス応答が短く、$N$ が小さかったので、$y_-(n)$ の畳み込みに FIR を使用したが、インパルス応答が長い場合は FFT を用いる方がよいだろう。
なお、理論的見通しをよくするため、本文では $y(n)$ の振幅スペクトルが全周波数で1になるように記述したが、$N$ が大きいほど $y(n)$ の振幅は小さくなるため、実際には測定系と被測定系の非線形性が問題にならない範囲で $y(n)$ に適当なゲイン $A\gt 1$ を掛けて大きくし、$y_-(n)$ または畳み込み結果にその逆数 $1/A$ を掛けて計算すべきである。
参考文献
- Borish, J. and Angell, J.B., "An Efficient Algorithm for Measuring the Impulse Response Using Pseudorandom Noise," J. Audio Eng. Soc., Vol. 31, No. 7/8, pp. 478–488, August 1983.
- 金田豊, 「M系列を用いたインパルス応答測定における誤差の実験的検討」, 日本音響学会誌, 52巻10号, pp. 752–759, 1996.
- 柏木濶, 『M系列とその応用』, 昭晃堂, p. 23.
- M系列 — Wikipedia
FAQ
M系列とは何ですか？
最大長系列（Maximum Length Sequence）の略で、擬似乱数系列の一種である。線形帰還シフトレジスタから生成され、自己相関がほぼ単位パルスに近い性質を持つ。この性質がインパルス応答測定に利用される。
なぜ時間反転した信号を畳み込むのですか？
振幅スペクトルが全周波数で1の信号とその時間反転信号の巡回畳み込みが単位パルス $\delta(n)$ になることを利用している。これにより、未知系の出力に反転信号を畳み込むだけでインパルス応答を復元できる。
先頭2周期が乱れるのはなぜですか？
測定開始前の入力が0であり、周期信号の周期性を前提とした計算と整合しないためである。$N$ サンプル同士の畳み込みが $2N-1$ サンプルになることから、先頭2周期では巡回畳み込みが成立しない。アベレージング時は先頭2周期（$2N$ サンプル）を除外する必要がある。