スプライン関数を使った補間について

現在位置のナビ

トップコンピュータの国雑記帳趣味趣味音響 → スプライン関数を使った補間について

説明

トラ技2018年10月号の45ページ第一章には、スプライン関数を用いた補間方式でPCMをオーバーサンプリングする記事が載っています。 筆者の溢れる情熱に対して、展開している技術論のレベルが追いついていません。 Webmasterがトラ技編集長に指摘したら、その時点ですでに他の読者から多数の指摘が来ていたそうです。 トラ技編集部には今更伝えるまでもありませんが、すでに届いた多数の指摘を見ていない他の読者のためにwebmasterの視点で指摘をします。

指摘の前に基本的な話

PCMをD/A変換する直前に、アナログLPFの設計を楽にするためにオーバーサンプリングする手法はオーソドックスです。 気をつけなければならないことは、オーバーサンプリングによって元のアナログデータが復活するわけでは無いことです。 デジタルデータに変換したことによって、PCMに格納されたデータは制限されています。 A/D変換の時に失われたデータを精密に復活する手法は存在しません。 今の所「理想的」と言える手法もありません。 一長一短です。

勘違いレベルの指摘

プレエコー

筆者はFIRのことをディジタル・フィルタと呼んでいて、プレエコーは、ディジタル・フィルタの欠点だと主張しています。 スプライン関数を用いても、デジタルデータの補間は一種のディジタル・フィルタです。 またプレエコーは、「補間の時に未来のデータを使用する」ときに出る特徴です。 このスプライン方式も補間に未来のデータを使用するので、プレエコーが出ています。 その証拠に、P.54図11(a)やP.45タイトルバックのグラフでも、sin波が始まる直前にグラフが一旦沈み込んでいます。

理想的な過渡応答波形

基本的な話で述べたように、どんな補間も理想的ではありません。 図1(a)で方形波が立ち上がる前の波型も立ち上がった後のリンギングもFIRのエコーとは限りません。 方形波は、周波数無限大の高域波形を含みます。 LPFで高域をカットしたら、立ち上がりは垂直から斜めに変わり、立ち上がり前後にリンギングが発生します。 方形波の元の波形を維持できないのは、高域が落ちているからであり、ポストエコーがどの程度影響しているかは、高域を落とした結果と分けて考えなければなりません。

SSDAC方式で方形波の波形が維持できるとすれば、ポストエコー以前の問題として、高域信号を十分に減衰できていない欠点を持つことになります。 過渡特性が良好なのではなく、LPFとして能力不足なのです。

PCM1704

記事中で、以下の記述があります。 「PCM1704(中略)設計も古く、新規設計には不向きです。」 PCM1704が新規設計に不向きなのは同意しますが、その理由は廃品種だから入手性に問題があるためです。 古い設計のチップがいまだにプレミア価格で取引されているのは、設計がこなれていて問題が無いからです。 古くても、良いものは良いのです。

イマドキのオーディオマニアは、「PCM1704はDSDを再生できない」とか言いそうです。 Webmasterに言わせれば、DSDを再生したかったら、705.6kHzのPCMに変換すればよいのです。

ニュートン方程式

P.46表2でSSDACの動作原理に「ニュートン方程式」とあります。 「ニュートン方程式」という表現は初めてみました。 よく似た言葉に「ニュートンの運動方程式」がありますが、SSDACは力学を扱ってはいませんね。 おそらく「時刻を変数にした高次多項式」が正確な表現でしょう。

量子化ノイズ

なんか、いろいろと勘違いしています。 「量子化ノイズ」とは、A/D変換するとき、デジタルで表現できる最小ステップからこぼれ落ちることになった微小量です。 オーバーサンプリングしても、元データの量子化ノイズを改善することはできません。 元データがCDDAの16bitだったら、補間の演算精度が24bitでも32bitでも16bit分の量子化ノイズを含みます。 P.47灰色の網掛け部分にある式 Nq = (Eq) 2 / 12をWikipediaの「量子化誤差」の記述と読み比べると、筆者がいくつも勘違いしていることがわかります。

無限スプライン関数を解いていない

記事中には「無限スプライン関数を解いた」とする記述があります。 AppendixのP.61式(21)以降を読めばわかるとおり、筆者は無限スプライン関数を解かずに、誤差が要求精度を下回った時点の有限回数で計算を打ち切っているだけです。

スプライン関数は元のサンプル点を必ず通る

筆者は、スプライン関数が元のサンプル点を必ず通ることを利点だと考えているようです。 ノイズシェーピングによってLSBに加減算が行われた場合は、LPFによってサンプル点を外れた波形が出てくるほうが嬉しいです。

有限スプラインの解き方

PCMデータは有限長です。 有限長であることを利用した解き方があります。 筆者は、「有限スプライン関数は、両端の点の2次微分値をゼロにすれば解くことができますが、」と言っています。 その方法では筆者自身が認めるように、全データを処理しなければなりません。 もっと簡単な解法があるので、後で書きます。

用語

筆者が何度も「1次微分」「2次微分」と書いているのは、「1階微分」「2階微分」の間違いのようです。 どうやら方程式の取扱には慣れていないようですね。 「均して」を「馴らして」と書き間違えるのはヒューマンエラーですから大目に見ますが、テクニカルタームを繰り返し間違えると、素人のように見えてしまいます。

アイソレータの場所

P.50右の段によれば、試作回路にはアイソレータが入っているそうです。 USB→ I2SのアダプタとFPGAの間にアイソレータを入れているようですが、アナログ回路を安定させる目的ならば、FPGAの後ろに入れたほうが良いのではないでしょうか。 出力クロックが高すぎて、適切なアイソレータが存在しなかったためだったりして。

ΔΣ

P.51左の段に、「D-Aコンバータは16ビットなので、残りの8ビットはΔΣまたは四捨五入で出しています」と書いています。 おそらく「ΔΣ」ではなくて「ノイズシェーピング」と書きたかったのでは無いでしょうか。

オーバーサンプリングをFIRで処理するとき元データに0を追加している

P.50図5(b)で、FIRフィルタに通す前のデータとして0データを追加しています。 この筆者に限らず、オーバーサンプリングの説明で同様のことを言う人が多いのですが、同じデータを並べて階段状にしたほうが、FIRのフィルタ係数精度を活かせるはずです。 事実、webmasterの作ったオーバーサンプリングフィルタではCDDAを8倍オーバーサンプリングする時に、同じデータを8個づつ並べています。 理論上も聴感上も、このほうが有利です。 計算量が増えるわけでもありませんし。 なぜ0データを追加する説明が広まっているのか不思議です。

FIRの欠点はそこじゃない

P.52のコラムで、「n(webmaster注 FIRのタップ数)を無限数にして、各点の値を図AのFIRフィルタのタップ係数に入れると、図Bに示すような完全な矩形の周波数特性のFIRフィルタになります」とあります。 筆者はwebmasterよりもFIRについての理解が足りないようです。 FIRはFFTを応用した技術であり、最大の欠点は「原理上周期的信号のみを対象としている」ことです。 それでは困るので、タップ数を調整したり、窓関数をかけたり工夫しています。 FIRも所詮近似のテクニックであり、100%完全な方法にはならないのです。

また、P.53の図DにはFIRの周波数特性なるグラフが載っています。 FIRの特性はタップ数と遮断特性で大きく変わるので、定性的な説明で誤解を増やしてほしくありません。

過渡応答の収束

P.54から引用します。 「過渡応答において、従来のディジタル・フィルタを使ったD-Aコンバータは、現在のデータから1データずつ遠ざかるにつれて、1/2、1/3、1/4、1/5……1/nというレートで影響力が薄れます」 ここ、webmasterには何を言っているのかわかりません。 P.53図Cに筆者が書いているように、FIRの係数は波打っています。 1/2 1/3 1/4…が並んでいるわけではありません。

計算量

P.56のコラムから引用します。 「SSDACは、1サンプリングの間に、219回の積和演算をしています」 Webmasterが2016年6月に第13回1bit研究会でデモンストレーションした機材は、IntelCPUがSSE命令を使って、352.8kHzfsのデータに8191タップのFIRを演算しています。 ちょっとした工夫をしているので44.1khzfsの1サンプルあたり、8x1025回の積と和を計算しています。 SSDACの10倍以上の計算量ですね。

さらに引用します。 「私のCorei5のパソコンは、ベンチマークでの実測値が1880MFLOPSでしたので、数字上は1313MFLOPSであれば処理できそうに見えますが、おそらく遅延時間が大きいためリアルタイムの処理は困難でしょう」 表記が抽象的なので、言いたいことがはっきりとはわかりません。 1bit研究会でのデモンストレーションでは、44100Hzx2chx8倍x1025tapx2(積と和)で1446MFLOPSの演算をしていました。 2016年6月当時は3GHzのCorei3を使っていましたが、Intelのドキュメントを読んでSSE命令を最適化したら、3GHzのPentiumDや1.7GHzのATOMアーキテクチャーでも余裕で動きます。 Webmasterの自宅では、リアルタイム処理はこんなんです。

フォン・ノイマン・ボトルネック

P.56のコラムから引用します。 「しかし、プロセッサにはノイマンのボトルネックという基本的な構造欠陥があり、演算1に対して100の命令実行が伴います」 ここも意味不明です。 おそらく、フォン・ノイマンのボトルネックの話をしているのでしょう。 「フォン・ノイマン型プロセッサは、命令コードもデータもメインメモリから持ってこなくてはならないので、そこが律速段階になる」という話です。 律速段階の話であって、構造欠陥ではありません。 また、演算1に対して100の命令実行とは、何を言いたいのかまたしても不明です。 マイクロコード・アーキテクチャーとボトルネックは別の話ですよね。 フォン・ノイマンのボトルネックは確かに存在しますが、SSE命令を使ったFIRなんて、命令コードもFIR係数もキャッシュに載るくらい小さいので、気にする必要は無いと思います。

2018年10月8日追記 3次方程式でsin()を近似することの妥当性

この章は解析的に論じますが、妥当性をきっちり評価せず、曖昧で終わります。 本来sin()波の重ね合わせである音楽を、3次関数で近似することが妥当なのか、できる範囲で論じます。 理工学部で習う数学の8割方を、読者がマスターしていることが前提です。 ここだけ文体が「だ、である」調になるのは、論文を書くときの癖なのでご容赦ください。

まず、基本です。

項が無限次まで続くことの考察

条件1: sin(x) のテイラー展開において

0<=x<=1 ならば
前の項ほど大きく後ろの項は0に収束する
1<xならば
前の項ほど小さく後ろの項は∞に発散する2021年4月25日改定 N番目の項はN < xだと増加傾向で x < N だと減少傾向

音楽信号を取り扱うので、x<0について考えなくても良い

条件1より、有限次の項までに計算を限定して後ろを切り捨てれば、1以下のxに関する大部分を保存しながら1以上のxの大部分を捨てることができる。 x の大小は信号周波数に比例するので、高次の項を切り捨てる効果はゆるい遮断特性のLPFとみなすことができる。 ただし、遮断特性がゆるすぎて希望しない信号(切り捨てたい高域が残ったり、低域の信号が変形したり)がでることも予想される。

ちなみに 3次までの項で計算を打ち切るということは ω 3 = 3! を満たす ω がカットオフ周波数となる。 ω について計算すると、ω = 1.817120593 となる。 もとのスプライン関数は、1サンプリング周期を区間[ 0,1 ]に正規化していたので、 ω =1.817120593 とはサンプリング周波数の ω / (2 π) 倍を意味し、CDDA ならば 12.8kHz になる。 本方式の LPF 効果は、遮断特性が非常に緩やかであることに注意。 ナイキスト周波数(ω = π)近辺を遮断周波数にしたい場合は、(π n)/(n!) がほぼ1 となるような n を探してn次関数としてスプライン関数を定義すればよい。

2018年10月8日追記 有限項数でスプライン関数を求める

もともとスプライン関数は、平面上のn個の点を滑らかな曲線で結ぶためのものである。 時間変化する関数に適用を試みると、末端の処理を考えなければならない。 再帰的に無限回の計算を試みるのは、無茶である。 PCMは有限個のデータ列なので、上手に扱えば末端から一意にスプライン関数が決まる。

PCM のデータ列 D 0, D 1, D 2, D 3, …が与えられたとする。 スプライン関数は3次方程式で、左右の区間との接続において関数値、1階微分値、2階微分値が一致するものとする。

方式考案者の計算方式は、区間 j=n を求めるために2区間 j=n-1, j=n+1のスプライン関数を必要とし、その前段階として3区間 j=n-2, j=n, j=n+2のスプライン関数を必要としていた。 再帰的であると同時に、前処理に次ぐ前処理を重ねていくと、無限個のスプライン関数を事前に求めておかなくてはならない。 漸化式と前処理の有限回数での打ち切りにより、計算量を限定していた。

新規に別の計算方式を提案する。 PCMのデータ列に新たに1個のデータ D -1 = 0 を追加する。 また、区間 j=-1 のスプライン関数を f -1( x ) = D 0 x 3 と定義する。 すると、区間 j=0 のスプライン関数は接続条件とD 1から一意に定まり、 f 0( x) = (D 1 - 10D 0 )x 3 + 6D 0 x 2 + 3D 0 x + D 0 となる。 j が1以上の区間も同様で、接続条件と次のサンプルデータ D {j+1} を用いて f j (x) がもとまってゆく。 もはや区間 j を求める際に、未来のデータ D {j+2}, D {j+3}, D {j+4}, … を必要としなくなる。

新規の計算方式は、スプライン関数の次数が4以上でも有効である。

御意見無様

このwebサイトは、読者の意見を歓迎します。 御意見無用ではなくて御意見無様(参考 下妻物語)です。 Webmasterが考え違いをしていたら、遠慮なくビシバシご指摘ください。 納得がいったら、ここで訂正します。

2020年2月11日追記 特許が査定されていない

トラ技の記述によれば、「本製品SSDACに使用されているデータ補間技術は、小林芳直により発明されSLDJ合同会社により現在特許出願中です。」とあります。 2020年2月4日現在、特許情報のデータベースで検索すると特許出願2017-130401 特開2019-016002が見つかりますが、査定種別は『査定なし』になっています。 出願したけど成立しそうもない特許を「出願したぜ」と威張っているのでしょうか。

2020年5月7日追記 審査請求出ました

J-PlatPatによれば、2020年4月16日にようやく審査請求が出たみたいです。 特許としては成立するかもしれませんが、Webmasterが提案した手法は特許に抵触せずにより簡単に使えますね。

2021年7月10日追記 拒絶理由通知出ました

J-PlatPatによれば、2021年7月6日に拒絶理由通知が出たみたいです。 「スプライン関数でオーバーサンプリングの補間を行うこと」に先進性はないという判断です。 「計算精度に応じて再帰計算を打ち切る」請求項には先進性が認められています。 さて、出願人から意見書は出るのでしょうか?

2021年7月19日追記 意見書出ました

J-PlatPatによれば、2021年7月16日に意見書が出たみたいです。

補正された請求項を読むと、「スプライン関数で補間しかつ、計算精度に応じて計算量を設定するシステム」として、新たに申請しています。 この特許が成立しても、webmasterがここで提案した計算方式はこの特許に抵触せず(他の特許に抵触するかどうかは自分で調べてね)、より少ない計算量で補間できます。

2021年7月31日追記 修整された特許が成立しました

J-PlatPatによれば、2021年7月30日に成立したみたいです。 特許料も8月19日に支払われています。 トラ技に載ったSSDAC方式をそのまま使用するには、特許使用料を支払わなければなりません。 Webmasterが提案した方式だったら、特許出願2017-130401には抵触しません。

更新日時

2018年10月7日 初出

2021年8月31日 追記


back button 趣味趣味音響へ