PythonWave

全体像

このページを一通り読んだら， waveファイル から1秒毎の \(L_{eq}\) をCSV形式で出力できるようになる。

おまけで，波形をグラフとして表示できるようになる。

まえがき

print('Hello World!!')

print(1 + 2)

準備

なにはともあれ，まずは環境構築。

Python

• ダウンロードページに進んで，"Windows x86-64 executable installer"をダウンロードする。
• 環境変数を通すようにチェックして，ダイアログにしたがってインストールする。

Visual Studio Code

• 作業フォルダを作る。
• Pythonのパスを通す。
• 拡張機能のインストール。

シンタックスハイライトからエラーチェック，

ライブラリが入ってないと警告してくれたり，とにかく便利。

2019/05/23時点最新版のPython3.7.3とVS Code1.34.0 zip版を使用する。

以下テストコード。これが動いたら今日から君もPythonist。

"好きなファイル名.py"にして保存，F5で実行。

calc.py

print("Hello World!!")
print(1 + 1) # 四則演算
print(2 - 1)
print(2 * 3)
print(4 / 3)
print(8 % 3)

import numpy # import: "このライブラリを使うよ"の宣言

arrFoo = numpy.asarray([1, 2, 3])
print(arrFoo)

import numpy as np # as: "この略で呼べるよ"の宣言

arrFoo = numpy.asarray([4, 5, 6])
arrBar = np.asarray([7, 8, 9])

print(arrFoo)
print(arrFoo[0])
print(arrFoo[1])
print(arrBar)

初回はnumpyがないよって怒られるはず。

Ctrl + @ でコンソールに移動。

pip install "インストールしたいライブラリ名"

でライブラリ追加。

pipが古いとか言われたら，指示に従ってインストール。

実現したいこと

A特性補正済みのwaveファイルから1秒毎の\(L_{eq}\)が知りたい。

分解する

( In ) waveファイル

• 構造理解

  音を記録するには？

  A-D Convert アナログ-デジタル変換

  PCM Palse Code Modulation パルス符号変調

  2つのパラメータ

  • sampling rate サンプリング周波数[Hz]

    1秒間に何回書き込むか。

  • quantization bit rate 量子化ビット数 (encoding bit rate 符号化ビット数)

    どのくらいの度数に分けるか？

  
  
• 取り込み方法

  Python nativeでいけるか，外部ライブラリはどれが使いやすそうか。

  標準ライブラリのwaveは拡張性がなさそう。

  モノラル/ステレオ以上の対応は書いてない．4chとかのアンビソニックス音源は読めない？

  読み出しても格納が面倒そう．24bit96kHzとかややこしげ。

  ライブラリはPySoundFileが良さげ。

  簡単，これに尽きる。

PySoundFileの使い方，matplotlibでグラフ表示。

plotWave.py

import numpy as np
import matplotlib.pyplot as plt
import soundfile as sf

# 各自書き換える
path = 'D:/documents/labo/pythonLeq/16bit1kHzMono.wav'

info = sf.SoundFile(path) # ファイル情報の確認
print('samplerate: ', info.samplerate)
print('channels: ', info.channels)
print('subtype: ', info.subtype)
print('format: ', info.format)

# データの読み込み
data, samplerate = sf.read(path) # data: 振幅(-1:1), samplerate: サンプリング周波数[Hz]
print(samplerate)
time = np.arange(0, len(data)) / samplerate # time: ファイルの長さ[s]
print('time: ', time)
print(time[-1]) # array[-1]: 最後の要素
fractionTime = len(data) % samplerate # 1秒未満の長さ
print('fractionTime: ', fractionTime)

plt.plot(time, data, label='wave')
plt.xlim(0, 0.003) # 表示する範囲
plt.ylim(-1, 1)
plt.xlabel('time(s)') # 軸ラベル
plt.ylabel('amplitude')

plt.show()
plt.close()

SPL sound plessure level, 音圧レベル[dB]を計算する

\(p_0\)は基準音圧[Pa]， (20 µPaのとき， dB SPL という。)

\(p_1\)は瞬時値[Pa]，

\(\log\) は \(\log_{10}\) : 常用対数を表す。

「ただし\(0\lt p_0\)，\(0\lt p_1\)」あるいは絶対値をとる。

calcSPL.py

import math

p0 = 20 * 10 ** -6
p1 = 20 * 10 ** -6 * 50000

print( 10 * math.log10(p1 ** 2 / p0 ** 2) ) # 数式と同じ意味
print( 10 * math.log10( (p1 / p0) ** 2 ) )
print( 20 * math.log10(p1 / p0) )

def calc(p0, p1): # 繰り返し使えるように関数宣言 define
  print(20 * math.log10(p1 / p0), 'dB')

print('twice')
calc(p0, 20 * 10 ** -6 * 2) # 2倍は？

print('10 times')
calc(p0, 20 * 10 ** -6 * 10) # 10倍は？

for i in range(1, 11): # 1から10倍
  print(i, 'time(s)')
  calc(p0, 20 * 10 ** -6 * i)

p0 = 20 * 10 ** -6 # = 0.00002 = 20 * 10 ^ -6

print('decimal point')
calc(p0, 20)
calc(p0, 2)
calc(p0, .2)
calc(p0, .02)
calc(p0, .002)
calc(p0, .0002)
print('bug?')
calc(p0, .00002)

print('use "for statement"')

for i in range(7):
  calc(p0, 20 * 10 ** -i)

print('test:', 20 * 10 ** -5)
print('test:', 20 * 10 ** -6) # != .00002

p0 = .00002
calc(p0, .00002)

今回の問題は"16bitが何dBか"ということ。

グラフを交えて実際に計算してみる。

graph.py

import math
import matplotlib.pyplot as plt

quantity = 2 ** 16
print((20 * math.log10(quantity / 1)))

x = []
y = []
for i in range(1, quantity):
  x.append(i)
  y.append(20 * math.log10(i / 1))
plt.plot(x, y)
# plt.xscale('log')
plt.grid(which='major', color='black', linestyle='-')
plt.xlabel('quantization bit rate')
plt.ylabel('SPL(dB)')
plt.show()
plt.close()

16bitが96dBとすると，これを超える音を録音できなくなる。

なので，\(p_0\)を"通常 20 μPa"ではなく可変にして，dB SPLに合わせる。

このとき， Calibration Signal 校正信号 から\(p_0\)を求める。

校正信号を基準にdB SPLを計算する。

dBを入力してもらって \(RMS\) から，

\(p_0=EU\): Engineering Unitとする。

leqSimu.py

import math

pCal = 0.5 # 校正信号のRMS値(計算結果)
lCal = 94 # 想定する校正信号の強度
eu = pCal / 10 ** (lCal / 20)
print(eu)
wavRms = 0.1 # 入力信号のRMS値(計算結果)
leq = 20 * math.log10( wavRms / eu )
print(leq)

import math
print(math.log2(10 ** (194 / 20))) # 音の限界

( Fn ) \(L_{eq}\)

音圧レベルの時間当たりの平均値。

わかるひとには，交流の利得の音バージョンっていえばわかる。

• 計算方法 \[L_{eq}=10\log\left[\frac{1}{t_1-t_0} \int_{t_0}^{t_1}\frac{{p_1}^2}{{p_0}^2}dt\right]\]

  \(p_0\)	: 基準音圧レベル
  \(p_1\)	: 収集した音圧
  \(t_0\)	: 測定の開始時間
  \(t_1\)	: 測定の終了時間

  参考

  Leq音圧レベルとは何ですか? NATIONAL INSTRUMENTS

  
  さっきの音圧レベル[dB]を，

  \[ L_{eq}=\class{mathbg}{10\log}\left[\frac{1}{t_1-t_0} \int_{t_0}^{t_1}\class{mathbg}{\frac{{p_1}^2}{{p_0}^2}}dt\right] \]

  その時間分全部足してあげて，

  \[ L_{eq}=10\log\left[\frac{1}{t_1-t_0} \class{mathbg}{\int_{t_0}^{t_1}}\frac{{p_1}^2}{{p_0}^2}\class{mathbg}{dt}\right] \]

  計測時間で割る。

  \[ L_{eq}=10\log\left[\class{mathbg}{\frac{1}{t_1-t_0}} \int_{t_0}^{t_1}\frac{{p_1}^2}{{p_0}^2}dt\right] \]
  

  \(RMS\): Root Mean Square value 二乗平均平方根

  = effective value 実効値

  \[RMS[x]=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(x_i)^2}\]

  testRMS.py

  import math
  
  def rms(array):
    ans = 0
    for i in range(len(array)):
      ans += array[i] ** 2
    ans = 1 / len(array) * ans
    ans = math.sqrt(ans)
    return ans
  
  x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
  print('calc:', rms(x))
  
  y = [1, 2, 3, -4, 5, 6, 7, 8, -9, 10]
  print('calc;', rms(y)) # マイナスがついてもok
  

  ファイルを取り込んで計算してみる。

  amplitude.py

  import numpy as np
  import soundfile as sf
  
  path = 'D:/documents/labo/pythonLeq/16bit1kHzMono.wav'
  
  data, samplerate = sf.read(path)
  time = np.arange(0, len(data)) / samplerate
  
  # RMSの計算
  # sqrt: 平方根
  # mean: 平均
  # [a: b]: aからbまで，ここでは1秒間
  rms1st = np.sqrt(np.mean(data[0: info.samplerate] ** 2))
  rms2nd = np.sqrt(np.mean(data[info.samplerate: info.samplerate * 2] ** 2))
  rms3rd = np.sqrt(np.mean(data[info.samplerate * 2: info.samplerate * 3] ** 2))
  print(rms1st, rms2nd, rms3rd)
  

  少し癖のあるnumpyの範囲指定

  setRange.py

  import math
  import numpy as np
  
  x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
  # print(x[0:(len(x) - 1)] ** 2) これはエラー
  
  y = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
  print(y, y ** 2)
  z = np.array(x)
  print(z, z[3:7] ** 2)
  print(z, z[0:9] ** 2)
  print(z, z[1:10] ** 2)
  print(z, z[0:10] ** 2)
  print(z, z[0:-1] ** 2) # 配列最後尾指定の常套句 -1 は使えない
  print(z, z[0:len(z)] ** 2) # 0番目から len(z) - 1 番目まで
  

  使えないnative mean

  uselessMean.py

  import math
  from statistics import mean
  import numpy as np
  
  x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
  print(math.sqrt(mean(x[0:len(x)] ** 2)))
  print(np.sqrt(np.mean(x[0:len(x)] ** 2)))
  
  print(sum(x[0:len(x)] ** 2))
  print(mean(x[0:len(x)] ** 2)) # 見事に小数点切り捨て
  print(np.mean(x[0:len(x)] ** 2))
  

  
  

  参考

  「収録した騒音計やマイクロホンからの音データをオースコープベーシックを使って音圧レベルとして表示する方法」 小野測器 - DR-7100 FAQ

( Out ) CSV

行列形式のテキストファイル。

今回の出力でお世話になる。

outputCsv.py

import datetime
import math
import numpy as np
testArray = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
# csv出力
filename = 'leq_' + datetime.datetime.today().strftime("%Y%m%dT%H%M%S") + '.csv'
np.savetxt(filename, testArray, delimiter=',', fmt='%d')

組み立てる

複数引数を取るためにArgumentParserを使用。

leq.py

import argparse
import datetime
import math
import numpy as np
import soundfile as sf

# 引数の設定
parser = argparse.ArgumentParser(description='Process some integers.')
parser.add_argument('inputFile')
parser.add_argument('calibrationFile')
parser.add_argument('SPL', type=float)
args = parser.parse_args()

# ( In ) データの読み込み
data, samplerate = sf.read(args.inputFile)
calData, calSR = sf.read(args.calibrationFile)
time = np.array(np.arange(0, (len(data) / samplerate)), dtype='int32')
print('time: ', time[-1])

# ( Fn ) 計算
def rms(amp, sec):
  return np.sqrt(np.mean(amp[samplerate*sec:samplerate*(sec+1)] ** 2))
rmsCal = rms(calData, 1)
leq = [20*math.log10(rms(data, i)/(rmsCal/10**(args.SPL/20))) for i in range(time[-1])]

# ( Out ) csv出力
filename = 'leq_' + datetime.datetime.today().strftime("%Y%m%dT%H%M%S") + '.csv'
np.savetxt(filename, leq, delimiter=',', fmt='%d')

あとがき

TODO

jupyter notebookとかgoogle colaboratoryにコードを移す

構成を考える。amplitude -> EU: rms, spl

csvはヘッダーと秒数くらいつける

複数ch対応

note Lmaxはnumpyのarray.max()とかですぐ出せる。