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