歪み率・周波数特性測定向けオーディオ信号発生の自動化（Python-Sounddevice）
2020/07/19, 2020/07/23, 2020/08/002 ページ作成中（URLも変更可能性あり）
WaveGene と WaveSpectra による歪み率・周波数特性測定で、WaveGene の「FFT用に最適化」の操作が非常に面倒くさいので、そういう周波数を WaveGene の代わりに自動進段で発生させることを試みました。手持ちの US-122MKii と WaveGene の相性が今一つよくない、という理由もあります。
周波数自動進段の間隔は広めにしています（変更可能）。学生実験で教わるように、周波数スイープをするときに最初からすべての周波数を総なめにするのではなく、まずは要所要所の特性を抑える、という考え方です。共振周波数近傍など細かいところを測定したかったら、あとで細かくやればよい。
こういうとき、零からプログラミングするのは今の時代、お馬鹿さん。Python にだいたい先人の例題があるので、ありがたく参考にさせていただきました。python-sounddevice というもので、あまりにも簡単にサウンドデバイスを制御できるとは、良い時代になったものです。
参考にしたサイト
- python-sounddevice 一般論
- ほぼ等比間隔の周波数設定
目次
- 初期検討
- 周波数振幅特性と歪み率特性の自動同時統計出力（Sansui STシリーズのひずみ測定）
初期検討
初期検討ソースコード
import sounddevice as sd
import numpy as np
# 使用するデバイス sd.query_devices() で得られるリストから適切なものを選択する。
# PC の起動や USB Audio I/F の接続変更のたびに異なるかもしれない。
dev = [12,12]
# サンプリングレート
fs = 96000
# WaveSpectra の FFT サンプルデータ点数
ps = 8192
# sin波の振幅
a = 1
# 発音するフレーム数（ps のこれ倍のサンプル数の波形を再生）
nframe = 8
# 上記のうちフェードインとフェードアウトに各々使うフレーム数
nframe_finfout = 2
# start ≦ n < stopで間隔はstep、(周波数分解能×これの要素値)Hzが出力周波数
mrange = np.arange(1, ps / 2, 64)
## ほぼ等比間隔の周波数設定にしたい場合：一旦 set に入れて重複要素を落とす
## mrange = list(sorted(set(np.geomspace(1, ps / 2, 64, dtype=int))))
# 現在の取得可能なデバイスの一覧を表示
print("query devices")
print(sd.query_devices())
print("周波数分解能", fs / ps)
# sounddevice の初期設定
sd.default.samplerate = fs
sd.default.channels = 2
sd.default.device = dev
# フェードインとフェードアウトのためのベクター
a_fadein = np.linspace(0.0, 1.0, num=int(ps * nframe_finfout))
a_fadeout = np.linspace(1.0, 0.0, num=int(ps * nframe_finfout))
for m in mrange:
# sin波の周波数
f = fs / ps * m
print("ステップ", m, ", 現在の周波数", f)
# sin波の wave table 生成
n = np.arange(ps * nframe)
swav = a * np.sin(2.0 * np.pi * m * (n / ps))
swav[0:int(ps * nframe_finfout)] *= a_fadein
swav[-1-int(ps * nframe_finfout):-1] *= a_fadeout
swav[ps * nframe -1] = 0
swav = np.append(swav, np.zeros(int(2*ps)))
# 再生
sd.play(np.array(swav), blocking=True) # 再生
sd.wait() # play()の処理が終わるまで待機
exit()
|
np.append で末尾にゼロを追加しているのは、短すぎるデータ列を再生しようとすると再生位置が正しく管理できずに途中で打ち切られてしまう場合があることへの暫定措置である。
特定の周波数帯を細かく測定するには、たとえば mrange = np.arange(1000, 1100, 2) などとすればよい。
サンプリングレートや WaveSpectra の FFT サンプルデータ点数は、使用するデバイスや測定目的に応じて適宜調整する。
発音フレーム数（上記の例ではフェードイン・フェードアウトを除いた振幅一定区間が 4 フレームになる）やフェードイン・フェードアウト長も、超低周波ではもう少し伸ばす必要があるかもしれない…が、そうやって見掛け倒しのデータをとってもしょうがない気もしないでもない。
初期使用
WaveSpectra の設定。US-122MKii は、なぜか、ASIO モードでは 32 bit を指定しないとうまく動かないことがある…。
US-122MKii の Line Out を Line In にループバックして実行中に、タイミングを選んでスクリーンショットを撮ったもの。矩形窓でもサイドローブが無いことを確認できる。WaveSpectra の FFT に対してちゃんと最適化されている。
出力の on/off に伴う過渡現象を軽減するため、先頭と末尾にフェードイン・フェードアウトを掛けている。それでも一般論として、過渡現象を含む区間を FFT すると数学的な性質で必ずサイドローブが現れるのは仕方ない。この現象はフェードイン・フェードアウトの時間を長くするほど軽減できるが、適当な匙加減で妥協する。
np.arange による等差数列の代わりに mrange = list(sorted(set(np.geomspace(1, ps / 2, 64, dtype=int)))) を用いてほぼ等比間隔の周波数設定にした場合。この設定では 117.1875 Hz 以下で周波数ステップが密接するため、ピークホールド線が横一線になっている。
これで OVERLAY.WSO ファイルを保存して、適当に加工すれば OK。WSOファイルは 8 要素ごとに改行されているので、上記の np.arange(1, ps / 2, 64) による等差数列の例では第 2 カラムの値を 8 行ごとに（8×8=64 要素ごとに）取り出せば、基本波応答振幅の周波数特性が描ける。
ほぼ等比間隔の周波数設定にした場合の整理には、ちょっと工夫がいる（未検討）。
周波数振幅特性と歪み率特性の自動同時統計出力（Sansui STシリーズのひずみ測定）
sd.play を sd.playrec にすれば同時録音もできる。また、今どんな周波数を出しているかはこのプログラム自身が知っている。よって、NumPy/SciPy の fft パッケージを併用すれば、その気になれば THD も RMS も出せるし、周波数振幅特性と歪み率特性の自動同時統計出力のようなことも、プログラミングの努力次第でできることになる。
…どうやらこれでいいらしい。基本波と2次〜9次高調波の振幅を出してみた。THD はいまのところ Excel で後処理計算することにする。
import sounddevice as sd
import numpy as np
import matplotlib.pyplot as plt
# 使用するデバイス sd.query_devices() で得られるリストから適切なものを選択する。
# PC の起動や USB Audio I/F の接続変更のたびに異なるかもしれない。
dev = [12,12]
# サンプリングレート
# Windows ではコントロールパネルのサウンドデバイスのレートと同じにしないとエラーになる
# （ASIO の場合はエラー後に再度実行しなおすとレート設定が反映されているらしい？）
fs = 96000
# WaveSpectraのFFTサンプルデータ点数、周波数分解能は fs / ps
ps = 8192
# sin波の振幅
a = 1
# 区分数
steps = 64
# 発音するフレーム数（ps のこれ倍のサンプル数の波形を再生）
# 切り出しの都合上、この値から 2*nframe_finfout を減じたものが最低2必要
nframe = 4
# 上記のうちフェードインとフェードアウトに各々使うフレーム数
# 切り出しを行うので0でもよいが、一応
nframe_finfout = 0.5
## # start ≦ n < stopで間隔はstep、(周波数分解能×これの要素値)Hzが出力周波数
## mrange = np.arange(1, ps / 2, steps)
# ほぼ等比間隔の周波数設定にしたい場合：一旦 set に入れて重複要素を落とす
mrange = list(sorted(set(np.geomspace(1, ps / 2, steps, dtype=int))))
# 高調波の解析上限次数
maxharmonics = 9
# 現在の取得可能なデバイスの一覧を表示
print("query devices")
print(sd.query_devices())
print("周波数分解能", fs / ps)
# sounddeviceの初期設定
sd.default.samplerate = fs
sd.default.channels = 2
sd.default.device = dev
# フェードインとフェードアウトのためのベクター
v1 = np.linspace(0.0, 1.0, num=int(ps * nframe_finfout))
a_fadein = np.sin(v1 * np.pi / 2)
#a_fadein = np.sqrt(v1)
a_fadein = v1
a_fadeout = a_fadein[::-1]
# 出力ファイルオープン
# ファイル名は決め打ち。自由に設定したければ sys.argv とかを使うように改良する。
file = open('test.csv', mode='w')
maxharmonics
file.write('freq')
for k in range(1, maxharmonics+1):
file.write(f", A{k}")
file.write('\n')
for m in mrange:
# sin波の周波数
f = fs / ps * m
print("ステップ", m, ", 現在の周波数", f)
# sin波の wave table 生成
n = np.arange(ps * nframe)
swav = a * np.sin(2.0 * np.pi * m * (n / ps))
swav[0:int(ps * nframe_finfout)] *= a_fadein
swav[-1-int(ps * nframe_finfout):-1] *= a_fadeout
swav[ps * nframe -1] = 0
swav = np.append(swav, np.zeros(int(2*ps)))
# 不平衡出力の Audio I/F で疑似平衡出力にするため Rch を Lch の逆相とする
# stack したものの転置（行がサンプル点数、列がch）を sd.playrec の再生入力とする
swav2 = np.stack([swav, -swav]).T
# 再生＆同時録音
record = sd.playrec(np.array(swav2))
sd.wait() # play()の処理が終わるまで待機
# 再生位置基準で一定振幅区間末尾の1フレーム分を切り出す
# sd.playrec では再生に対して録音が若干遅延するので、これでok
record1 = record[int(ps*(nframe-nframe_finfout-1)):int(ps*(nframe-nframe_finfout))]
#print(len(record))
#print(record)
#print(record1)
#print(record1.T)
# L,R に分離
record1L = record1.T[0]
record1R = record1.T[1]
# Lch のみ FFT 処理、Rch はとりあえず放置
print('fft... ', end="")
#print(len(record1L))
#print(record1L)
recfftL = np.abs(np.fft.rfft(record1L)) / ps * 2
print('fundamental_freq=', f, ', A1=', recfftL[int(m)], 20*np.log10(recfftL[int(m)]) )
value = str(f)
# 基本波と2次〜9次高調波の振幅を可能な範囲で出力
for k in range(1, maxharmonics+1):
if m*k < ps/2:
print(k, ', ', end="")
value += ',' + str(recfftL[int(m*k)])
print()
value += '\n'
file.write(value)
#plt.subplot(3, 1, 1)
#plt.plot(swav2, label="played sin wave") # 再生した信号
#plt.legend()
#plt.subplot(3, 1, 2)
#plt.plot(record1, label="recorded sin wave") # 録音された信号
#plt.legend()
#plt.subplot(3, 1, 3)
#plt.plot(recfftL, label="FFT Lch") # 録音された信号の Lch の FFT
#plt.legend()
#plt.show()
exit()
|
ASIO モードとして定義されているデバイスを選ぶのがみそで、MME はソフトクリッピング特性を持たせてあるのでこの目的には不適、WASAPI もいまひとつ。また、Windows のシステムサウンドが出るとサウンドデバイスの設定がリセットされてしまうようだ。その場合、このスクリプトをエラーに関わらずしつこく実行すると通る場合もある。
このスクリプトを開発したのは、ぺるけ式バランス型ヘッドホンアンプに内蔵する AKI.DAC のトランス探しが目的。
US-122mkii の Line Out Lch を Line In Lch に直結ループバックを 650mV で行った。INPUT LEVEL つまみを調整して -0.3dB F.S. くらいになるようにしている。Line In に不平衡入力を与えているためか、二次歪みが多い。Line Out Rch にせっかく出している逆相出力をちゃんと活用すれば改善されるのかもしれないが、まだできていない。また出力レベル、INPUT LEVEL の調整次第で傾向が結構変わる。そんなに高級な測定はできないが、割り切って使うことにする。
US-122mkii の出力を、そのまま直結と NJM4556A 反転1倍したもので平衡ループバックとしたもの。振幅を US0122mkii の PHONES/LINE OUT で絞り、片側あたり 325mV で Hot-Cold 間トータル 650mV。2次ひずみが大きく減り、3次ひずみが主体である。一方、数 kHz 以上では NJM4556A のひずみが見え始めているのかもしれない。
US-122mkii の出力を NJM4556A 反転1倍で低インピーダンス化して、LHLC08NB 1mH×2（磁気結合があるように配置したので2.3mH位あると思う）、13600pF、330Ωにてフィルタをかけ、山水ST-27の二次側（2kΩ側）に650mVとした場合の、20kΩ側出力の各次数レベルと、歪み率（レベルを基本波で規格化して100倍した数値…THD-Fに相当、ノイズは除外）。トランス二次負荷は US-122mkii ライン入力の Zi=15kΩ である。平衡入力にすると US-1222mkii 自身の二次歪みは気にならなくなる。11.71875 Hz で 1% 程度、その一つ上の 23.4375 Hz で 0.396% なので、決して優秀ではないが、実用に堪えないというほどでもないだろう。スクリプトで測定周波数の等比数列的間引きを行っているので、この測定が（回路の準備を除くと）2分ほどで出来てしまう。全部の周波数をリニアに愚直にやるのに比べて大幅に効率が良い。高調波が積算できる範囲は周波数が高くなるほど狭くなり、15kHz付近以上では3次高調波が取れないため見かけ上のDall(THD)が大きく下がっているが、本当に低歪みなわけではないので注意。
ぺるけ式バランス型ヘッドホンアンプでは、2V F.S. の I/F 出力を受けると自分の使い方ではボリュームが9時でも音量が大きいくらいに感じるので、昇圧無しの等倍で十分なようだ。そこで同上で山水ST-73Aの一次側（センタータップがない側）に650mVとした場合。INPUT LEVEL の増幅度はこちらのほうが大きいため、歪み率やノイズはわずかに不利になる。11.71875 Hz で 2.1% 程度、その一つ上の 23.4375 Hz で 0.57% なので、超低域は少し厳しい。しかし 100Hz 程度から上は ST-27 とさほど変わらない。400Hz 前後がこちらのほうが素直なのは US-122mkii の入力初段がバイポーラトランジスタ (2SA1084) であることも理由かもしれない。まあ一長一短といったところ。
適合インピーダンス表記が違うが巻線比が同じ ST-23 や ST-78 でも試してみたいところ。ステップダウンになる ST-29A、ST-25A、ST-76A なんかも検討の余地があるかもしれない。
なお、ゲルマニウム時代からある20〜40番台といくぶん時代が下った70〜80番台との間には、設計思想の相違があるのではないかとなんとなく思っている。
LCR ローパスフィルタ無しにて、山水ST-27の本来の信号通過方向 20kΩ側に 650mV および 325mV 印加したときの 2kΩ 側出力の観察。2kΩ 側センタータップは GND に接続。325mV 時は、振幅を US0122mkii の PHONES/LINE OUT で絞っている。LINE IN フェーダーの設定は同じ。
650mV では US-122mkii の LINE IN での増幅の歪みがかなり入っていて、ST-27 の本来の実力ではなさそう。（※US-122mkii フロントエンドは 2SA1084 エミッタ間の抵抗でゲインを決める回路のようで、高ゲイン状態で入力レベルが大きいと歪みが急激に大きくなる。）
325mV がほぼこの使い方の実力だろう。しかし、飽和しないはずなのに、歪み率がなかなか下がらない。5次や7次が長く尾を引いているのも気になる。これらは、コア自身の歪み、ならびに、巻線抵抗が増えることでコアの非線形性を打ち消すような十分な励磁電流（非線形の透磁率に対して励磁電流波形が歪むことによって、鉄心内の磁束の磁束の歪みが減る）を流せなくなることからくる歪みであろう。山水STシリーズの CB-19 系列は 12Hz にて 1% あたりが限界で、5次や7次とのバランスでは 1kΩ〜3kΩ 程度のものが AKI.DAC の出力レベルにマッチするようである。ただし、一部にインピーダンス公称値に対して巻線直流抵抗が妙に低い機種（具体的には ST-24 の 1kΩ 側、 ST-75 の 600 Ω側、 ST-76A）があり、なんともいえない。
略歴
- 2020/07/19 思いつきから3時間くらいでこのページを作り公開。
- 2020/07/23 先頭と末尾にフェードイン・フェードアウトを掛けるよう改良。ほぼ等比間隔の周波数設定の試行。
- 2020/07/25 基本波と2次〜5次高調波の振幅をファイルに書き出す試行。
- 2020/07/26 微修正、Lch-Rchでの疑似平衡出力対応。
- 2020/08/01 入力フルスケールに対する表示振幅を1に正規化、コンソールへのdB表示を追加。
- 2020/08/02 山水ST-27 と ST-73A の測定例。
- 2020/08/09 平衡ループバックと、山水ST-27 本来の信号通過方向での測定例。