最近のアンプは安くてもキャリブレーション機能があってすごい
こんにちわ、こんばんわ。アンプが壊れたのに買ってもらえない、かえるのクーの助手の「井戸中 聖」（いとなか セイ）でございます。
さて、いろいろな音響関連プログラムを作成すると、入力ファイルにホワイトノイズやピンクノイズを指定して、音響特性を検証することはよくあります。
たま～に、単一周波数(Sin波)で音響特性の計測をしたいときがありますが、その都度校正用の音ファイルを作成していました。簡単なので、すぐできるのですが、ナイキスト周波数ぎりぎりまでの波形作成は結構面倒（振幅が安定しないなど）なことがあり、Pythonプログラムで波形を作っておくことにしました。
波形の仕様
・16Bit 18Hz～24000Hz スィープsin波形/48K （とりあえず16Bitで実験することが多いため16Bit）：いわゆるDVD/衛星放送 音質です。
・長さ 上昇 2分、下降2分で計 4分
・振幅最大は Bit最大値 - 1
・開始、終了の 0.1 Sec は、フェードイン、フェードアウトする
※その他：24000Hz付近の振幅にはこだわらず、上昇、下降の反転ポイントで位相がつながるようにします
周波数をスィープさせながら、所定のポイントで正しい周波数にするのは結構難しいので、「校正用ファイル」自体の補正は実測で調整しました。
開始部分
0.1 Secの立ち上がりでフェードインしています。開始は20Hzくらいです
上昇・下降折り返し部分
ナイキスト折り返し部分で、位相の不連続はありません。120秒点で（ほぼ）24KHzとなるように調整しています。
終了部分
末尾0.1 Secでフェードアウトしています。終了も20Hzくらいです。
周波数特性（スペクトル）
20Hz部分でフェードイン/アウトしているので、20Hz付近は少し落ち込みます。
24KHz付近まで、きれいに -3dB/oct の直線を描いています。
振幅が同じであれば、スィープ周波数により -3dB/oct の特性となります（これで正しいです）
※オーディオインターフェースによっては（周波数/特性の相性やCPU負荷状況により）少しノイズが乗る場合があります。
実際の音
聞こえる音量感よりもパワーがありますので、聞く場合は音量は控えめ（通常の1/5～1/10程度）で聞いてください。周波数が高くなって聞こえなくなっても音は鳴っていますので、決して音量を上げないでください。
これで、いろいろな音響測定に使えます。（主にウェーブレットデコードの校正に使うつもりです）SoundCloudはアップロード時にリサンプリングされるようで、高い周波数は歪まくりでした。アップしたファイルは使用に耐えないと思いますので、もし使われたい場合は以下のプログラムを実行して、ファイルを作成くださいませ。（SoundCloudはあんまりなので44.1KHz版に差し替えるかもしれません）
プログラム
#! /usr/bin/python
# -*- coding: utf-8 -*-
import math
import numpy as np
import wave
CALIB_FILE = "calibration20to24K.wav"
SAMPLE_RATE: int = 48000 # サンプリングレート
WAVE_SEC: int = 240
AMP_MAX: int = 32767
FADE_INOUT: float = 0.1
def write_soundData(stream, filename):
s_write = wave.open(filename, 'w')
s_write.setparams((1, 2, SAMPLE_RATE, 0, 'NONE', 'Uncompressed'))
s_write.writeframes(stream)
s_write.close()
# WAVEフォーマットファイルの出力
def output_wave_file(x_hat, file_name, amax):
decode = x_hat / amax * AMP_MAX # 16bitでほどよく収まるように振幅を調整する
rtn_data = decode.copy()
decode.astype('int16') # Wave出力用に16Bit化します
work = b''
out_stream = b''
amp = 0
print('*** convert Array to Stream Now ***') # WAVEファイルとして出力します。
for i in range(len(decode)):
amp = int(decode[i])
if amp > AMP_MAX:
amp = AMP_MAX # クリッピング
if amp < -AMP_MAX:
amp = -AMP_MAX # クリッピング
ampbyte = amp.to_bytes(2, byteorder='little', signed=True) # 地道に1サンプルづつ変換します。
work = b''.join([work, ampbyte])
if i % 10000 == 0: # 一旦ワークに入れるのはそのままjoinしていくと指数的に遅くなるからです。
if i % 100000 == 0:
print('stream:%d' % i) # 遅いので、経過をレポートします。
out_stream = b''.join([out_stream, work])
work = b''
out_stream = b''.join([out_stream, work]) # 最後の残りをフラッシュします。
print('stream:%d' % i)
print('*** output Now ***')
write_soundData(out_stream, file_name) # ファイルに書き出します
print('DecodeFile = ' + file_name)
return rtn_data
def get_freq(asec, start, end, x):
eunit_log = math.log2(end / start)
ratio = 1.0 - abs(1.0 - 2 * (x / asec))
x_log = eunit_log * ratio
freq = start * 2 ** x_log
return freq
def genarate_wave(asec):
FREQ_ADJ = 0.441005 # 実測による調整パラメータ（折り返し点でナイキスト周波数限界付近になるように調整）
start_freq = 20
end_freq = SAMPLE_RATE / 2
x = np.linspace(0, asec, SAMPLE_RATE * asec)
freq = get_freq(asec, start_freq, end_freq, x)
adj = (x >= asec / 2) * (-1.0) + (x < asec / 2) * 1.0
fade = (x < FADE_INOUT) * (x * 10)
fade += (FADE_INOUT <= x) * (x<= (asec - FADE_INOUT)) * 1.0
fade += ((asec - FADE_INOUT) < x) * (asec - x) * 10
y = np.sin(freq * asec * FREQ_ADJ) / 2 * adj * fade
return y
if __name__ == '__main__':
wv = genarate_wave(WAVE_SEC)
amax = max(np.max(wv), -np.min(wv)) # 振幅の最大値を求める
output_wave_file(wv, CALIB_FILE, amax)
