この記事では、これまで紹介してきたSSDACを実際にパソコン上のPythonで実装し、デジタルデータをスプライン関数で補間するという作業をご理解いただきたいと思います。
今回は、まずは簡単な信号を録音したCDと同じ44.1kHz16bitのWAVファイルを読み込んで、任意のスーパーサンプリング倍率でスプライン補間した場合の波形をシミュレーションします。
最終的にはPythonで実装したソフトで、CDのデータをスプライン関数補間によってSRC（サンプルレートコンバート）したファイルを出力し、疑似ハイレゾファイルとして再生できるようにしたいと思います。
なおこの記事を書くにあたっては、スプライン関数を導出するための公式を使っています。スプライン関数の公式の詳細な説明は次の文献をご参照ください。
・トランジスタ技術2018年10月号記事
・電気学会論文 ECT-17-058-3
・電子情報通信学会論文 ICD2017-28
1．なぜスプライン関数で補間したいのか
オーディオDACでは通常、オーバーサンプリングという技術でデータ補間を行っています。これはごくおおざっぱに言うと、上限周波数20kHzに対してサンプリング周波数が44.1kHzでは粗いので、補間してなめらかにしてから再生するということです。また、このときオーバーサンプリングによって量子化ノイズやエイリアシングノイズの周波数を高域に追いやることで、ポスト・ローパスフィルタの設計が楽にできるというという大きな利点があります。
通常、オーバーサンプリングは、データに対して補間点に0値の点を挿入（インターポーレーション）して、FIRフィルタを通すことで実現していますが、ここで使われているFIRフィルタには、過渡的な再生信号に対してプリエコー、ポストエコーと呼ばれるノイズが発生するという、宿命的な欠点があります。図1に示すように、ポストエコー、プリエコーはサンプリング周波数の半分、たとえばCDだと44.1/2＝22.05kHzで発生します。
それでは、プリエコー、ポストエコーの発生を抑えつつ補間する方法はないか？ということで考案されたのが、スプライン関数による補間です。図2に示すように、スプライン関数による補間では、プリエコー、ポストエコーが大幅に抑えられているのがわかります。
これこそが、SSDAC最大の利点なのです。
図1．ESS9018K2Mによるポストエコー、プリエコー
図2．スプライン関数補間により、プリエコー、ポストエコーが抑えられる
2．なぜこれまでスプライン関数補間によるDACがなかったのか
飛び飛びの点をなめらかにつなぐために、たとえばEXCELでは包絡線の機能があり、これはスプライン関数によるものです。離散的な点をスプライン関数でつなぐ、という発想はごく当然なものです。それではなぜこれまでデジタルオーディオで応用されなかったのでしょうか。
スプライン関数を正確に計算するためには、開始点と終止点が確定していることが必要条件です。たとえば音楽再生で考えると、スプライン関数を計算するためには、ある曲が始まってから終了するまでのすべてのデータを使って演算する必要があるのです。そうすると、データを一曲分先読みしておいて、スプライン関数をすべて計算してから再生するということになり、再生ボタンを押してから再生開始まで、スプライン関数の演算時間を待たなければなりません。短い曲ならいいかもしれませんが、長い曲ほど演算に時間がかかります。また早送りなどをして途中から再生したい場合などには、全体を通して再生する場合とスプライン関数の演算結果が違ってしまうなどの不都合が生じます。
SSDACでは、これらの問題点について数学的なアプローチで解決することに成功しました。
3．どのように解決したのか
いま、CDからデジタル信号が入力しているとして、これは44.1kHz16bitのサンプリングなので1/44100≒22.68μｓごとに-32768～32767までのいずれかの値を持つデータが1個（ステレオだと2個）送られてきますから、これをアナログの値に変換するのがD/Aコンバータの仕事です。-32768～32767の値を持ったデータは、D/Aコンバータの仕様によって、たとえば0V～2Vとか、±2.5Vとかの範囲の電圧に変換されて出力されます。
今現在D/A変換している生データをD0として、ひとつ未来の生データをD+1、n個未来の生データをD+n、ひとつ過去の生データをD-1、n個過去の生データをD-nというように、現在の生データD0を中心に、データに番号を振っておきます。
スプライン関数でデータを補間するとは、具体的には現在の生データD0から1個未来の生データD+1までの間を、
・・・・・①
という3次関数でなめらかにつなぐということです。この式のd0は現在の生データD0です。
このときxの変域は
であり、式さえ求まってしまえば、あとはxを細かくステップすればするほどなめらかな補間ができます。
この式でｘ＝0のときはｙ0＝d0＝D0、つまり現在の生データ、x＝1のときは1個未来の生データD+1となります。
それでは、スプライン関数の公式を使って、式１を決定するやり方を見ていきましょう。式①を決定するとは、各係数a，ｂ，cを決定することです。
まず式①で、x^2の項の係数bを求める公式は、
・・・・・②
djは入力された生データです。e，fは次式で与えられます。
・・・・・③
・・・・・④
ここで、αは、
という定数です（これは発明者の小林芳直氏によってIS（Infinite Spline）定数と名付けられました）。
上記式③、④および定数αが、SSDACでスプライン関数の計算が可能になった最大のキモなので、以下に簡単に説明します。
スプライン関数の計算を実装可能にしたのは、このαという定数の発見と、それに伴い計算の打ち切り条件が算出されたことです。
これらは上述のトランジスタ技術2018年10月号にて詳しく解説していますので、ここではごくおおざっぱに説明します。
まず、現在の生データとひとつ未来の生データの間を補間するスプライン関数に与える前後の生データの影響は、1つ未来または1つ過去に遠ざかるにつれて、αを乗じた分だけ影響度が低下してゆく、ということが発見されました。つまり、式③、式④にて示されるように、ひとつ過去また未来の生データの影響度はα倍、2つ過去または未来の生データの影響度はα^2倍・・・ｎ個未来の過去または未来の生データの影響度はα^n倍となっています。αはおよそ-0.26ですから、2乗、3乗、4乗……と乗じるごとにどんどん小さくなっていきます。
さて、式③、式④では未来または過去に13個離れた生データまでで打ち切られています。これはどうしてかというと、これらの式はデータ精度を24bit出力として得る場合のものなのですが、13個より遠い生データが現在のスプライン関数の演算結果に与える影響は24bit精度以下になるので14個目以降の生データを打ち切ることができる、ということを示しています。
どこで打ち切っていいのか、ということは、希望する出力のbit精度によって違い、次のようになります。
16bit・・・・・未来過去それぞれ9個の生データまでで、10個以降は打ち切り
24bit・・・・・未来過去それぞれ13個の生データまでで、14個以降は打ち切り
32bit・・・・・未来過去それぞれ17個の生データまでで、18個以降は打ち切り
つまり、得たい出力ビット数に応じて、ビット数/2+1までの未来過去までを演算すれば良い、ということがわかったのです。
ここまででスプライン関数のx^2の係数bが計算できたので、残りのaとcを、次の公式を使って求めます。
・・・・・⑤
・・・・・⑥
ここで、djは現在の生データ、dj+1は1つ未来の生データ、bjは先に計算した現在のスプライン関数の係数b、bj+1は1つ未来のスプライン関数の係数bです。
これで、現在の生データから1つ未来の生データまでを補間するスプライン関数が求まりました。
4．PC上のPythonで実装する
それでは以上に述べたスプライン関数によるデータ補間を、実際にPythonで実装してみましょう。
ここでは、CDと同じ44.1kHz16bitサンプリングのWAVファイルのデータを、スプライン関数によって補間し、補間前後の波形をグラフで表示して比較する、ということを行います。
まず筆者の作業環境ですが、windows11にAnaconda3をインストールするとJupyter Notebookが自動的にインストールされるので、これを使っています。
Jupyter Notebookを起動すると、ブラウザ上にjupyterのHomeが起動するので、ここで作業を行います。画面右上のNewのプルダウンからNotebookを選ぶと、新しくNotebookが開きますので、ここのコード入力窓にコードを入力して最終行末尾にカーソルを持って行ってshift+Enterすれば実行されます。
Pythonを使い慣れた人はanacondaに付随しているSpyderを使う方が便利かもしれません。
次のコードの下にソースコードのダウンロードリンクを張っておきます。
以下に、今回Pythonで作成したSSDACのコードを示します。
##############ここから#################################
from os.path import dirname, join as pjoin #文字列接続メソッドjoinのインスタンス定義
from scipy.io import wavfile
import scipy.io
import numpy as np
import matplotlib.pyplot as plt
from scipy.io.wavfile import read
data_dir = 'C:\\MyData\\Python'
fname = '1k_5ms.wav'
wav_fname = pjoin(data_dir, fname)
alpha=-2+3**0.5 #IS定数
ss_ratio=4 #スーパーサンプリング比
samplerate, data = wavfile.read(filename=wav_fname)
#24bit精度でスプライン関数を計算するため、データの先頭と末尾に(24/2)+1=13個の計算用ダミーデータを追加する。
padding_array = np.array([[0, 0], [0, 0], [0, 0], [0, 0],[0, 0], [0, 0], [0, 0], [0, 0],[0, 0], [0, 0], [0, 0], [0, 0],[0,0]])
data=np.vstack((padding_array,data)) #data先頭にダミーデータ13個追加
data=np.vstack((data,padding_array)) #data末尾にダミーデータ13個追加
#ファイル全体分のインパルスe,fの空の配列（L,R）を用意する
e = np.zeros(shape=(len(data),2))
f = np.zeros(shape=(len(data),2))
#ファイル全体分のスプライン関数のaj,bj,cjの空の配列（L,R）を用意する
aj = np.zeros(shape=(len(data),2))
bj = np.zeros(shape=(len(data),2))
cj = np.zeros(shape=(len(data),2))
#仕上がりデータ
y_all=np.empty((1,2)) #ここで空の配列を作る
#ファイル全体分のインパルスe,fを計算する
for i in range(13,len(data)-13):
for j in range(0,14):
e[i]+=(alpha**j)*data[i-j]
f[i]+=(alpha**j)*data[i+j]
#先頭のe(-1)=0
e[12]=0
#末尾のf(+1)=0
f[len(data)-13]=0
#ファイル全体分のスプライン関数の係数aj,bj,cjを計算する
for i in range(13,len(data)-13):
bj[i]=-3*(3**0.5-1)*data[i]+3*(2*(3**0.5)-3)*(f[i+1]+e[i-1])
for i in range(13,len(data)-13):
aj[i]=(bj[i+1]-bj[i])/3
cj[i]=data[i+1]-data[i]-2*bj[i]/3-bj[i+1]/3
for i in range(13,len(data)-26):
for j in range(0,ss_ratio):
y_all=np.vstack((y_all,aj[i]*(j/ss_ratio)**3+bj[i]*(j/ss_ratio)**2+cj[i]*(j/ss_ratio)+data[i]))
#plt.plot(y_all)
x=np.arange(0,len(data),1)
plt.xlim(24,55)
#(24,55)@1k_5ms
#(24,35)@10k_500us
plt.step(x,data[x])
plt.title(fname+" NOS Zoom in")
plt.show()
x=np.arange(0,len(y_all),1)
plt.xlim(220,350)
#(110,170)@1k_5ms_ss2
#(220,350)@1k_5ms_ss4
#(900,1400)@1k_5ms_ss16
#(0,17)@10k_500us_ss2
#(0,34)@10k_500us_ss4
#(0,80)@10k_500us_ss16
plt.step(x,y_all[x])
plt.title(fname+" SS"+str(ss_ratio)+"x"+" Zoom in")
plt.show()
x=np.arange(0,len(y_all),1)
plt.xlim(0,len(y_all))
plt.step(x,y_all[x])
plt.title(fname+" SS"+str(ss_ratio)+"x"+" All over")
plt.show()
########################ここまで########################
ダウンロード - 20240716python_ssdac000.py
今回、コードで使用するための44.1kHz16bitサンプリングの短いWAVファイルを用意しました。これはefu氏によるWaveGeneによって生成したもので、Lチャンネルに方形波、Rチャンネルに正弦波を録音してあります。
ダウンロード - 1k_5ms.wav
ダウンロード - 10k_500us.wav
まず、コードの9，10行目で、WAVファイルのディレクトリとファイル名を設定しています。
data_dir = 'C:\\MyData\\Python'
fname = '1k_5ms.wav'
ここはご利用の環境に応じて適宜変更してください。
次に14行目でスーパーサンプリング比を設定しています。これは現在の生データから次の生データまでの間を何倍のデータで補間するかを設定します。現在頒布しているSSDAC基板では、SSDAC128なら128倍、PCM1704やI2Sデバイスを使ったSSDACは16倍になっています。
ss_ratio=4 #スーパーサンプリング比
66行目以降はグラフ表示の部分で、ここでは入力信号の周波数とスーパーサンプリング比によって適切なズーム倍率で表示されるように、plt.xlimの値を変更する必要があります。これは注釈としてコード中に記載してあるので、適宜変更します。
なお、グラフ表示にあたっては、実際のD/Aコンバータと同じく0次ホールドでグラフを描画しているため、ステップ表示になります。
その他の注意点としては、上で述べたとおり、24bit精度でスプライン関数を計算するために、現在のデータの前後13データが必要になるため、ファイルから読み込んだデータの先頭と末尾に13個の0データを付け加えています（18～22行目）。
それではさっそく実行してみましょう。
まずは1k_5ms.wavを使って、4倍スーパーサンプリングの波形を見てみます。ファイルのディレクトリだけ注意して、あとは上のコードのままでいけるはずです。次のグラフが表示されれば成功です。
図3．1kHzの方形波および正弦波のNOS（元）データ
図4．1kHzの方形波および正弦波の4倍スーパーサンプリングデータ（ズーム）
図5．1kHzの方形波および正弦波の4倍スーパーサンプリングデータ（全体）
図3が元のデータ、いわゆるNOSでのグラフで、図4が4倍スーパーサンプリング処理をしたものです。
図4ではデータが4倍スーパーサンプリングされて、ギザギザ（量子化ノイズ）が細かくなっている様子がわかります。
次に、同じ1kHzの信号を16倍でスーパーサンプリングした波形を見てみます。
プログラムは、14行目
ss_ratio=16 #スーパーサンプリング比
と75行目
plt.xlim(900,1400)
のように変更します。
図6．1kHzの方形波および正弦波のNOS（元）データ
図7．1kHzの方形波および正弦波の16倍スーパーサンプリングデータ（ズーム）
図8．1kHzの方形波および正弦波の16倍スーパーサンプリングデータ（全体）
図7を見ると、ほとんど段差が見られないほどなめらかになっています。
次に10kHzの信号ではどうでしょうか。10行目を
fname = '10k_500us.wav'
14行目を
ss_ratio=4 #スーパーサンプリング比
とし、67行目を
plt.xlim(24,35)
75行目を
plt.xlim(0,34)
に変更して実行します。
図9．10kHzの方形波および正弦波のNOS（元）データ
図10．10kHzの方形波および正弦波の4倍スーパーサンプリングデータ（ズーム）
図9の元データでは、正弦波の面影はなくなってしまっていますが、図10の4倍スーパーサンプリングではかろうじて正弦波が復元されているように見えます。ここでは全体データは図10とほぼ同じくなるため省略します。
次に同じ10kHzの信号を16倍スーパーサンプリングします。14行目を
ss_ratio=16 #スーパーサンプリング比
に変更し、75行目を
plt.xlim(0,80)
とします。
図11．10kHzの方形波および正弦波の16倍スーパーサンプリングデータ（ズーム）
元のデータは図9と同様です。
16倍スーパーサンプリングでは図10と比べてかなりなめらかになりました。全体波形は図11とほぼ同じなため省略します。
スプライン関数によって信号補間を行うSSDACの動作がご理解いただけたでしょうか。
このプログラムにより、今回ご用意した1kHz、10kHzの正弦波、方形波の他にも任意のWAVファイルを処理することができ、スーパーサンプリング倍率も任意に設定可能です。ぜひいろいろな波形、いろいろな条件で試していただきたいと思います。
さて、この記事の続編として、今回と同じくPythonで記述したソフトによってスーパーサンプリングした信号をファイルとして取り出し、擬似的にハイレゾ化する検証をする予定です。つまりスーパーサンプリングによるSRC（サンプリングレートコンバーター）です。
通常、アップサンプリングしたからといって音質が良くなるわけではなく、条件によってはFIRフィルタの使用によってプリエコー、ポストエコーが付加されるなど、百害あって一利無し、と考えられますが、スーパーサンプリングによってアップサンプリングし、それを市販のDACで再生した場合、発生するプリエコー、ポストエコーの周波数はスーパーサンプリングの倍率分高くなりますので、聴感上の悪影響が減少すると考えられます。
関連記事
SSDACをパソコン上のPythonで実装する ②
SSSRC（ スーパーサンプリング・サンプリングレートコンバータ ） ソフトウェア公開
最近のコメント