Blog

Python3でPink-TSP信号を生成

音響機器や室内音場の特性を調べる時には、ほとんどの場合対象のインパルス応答を測定します。

例えば自分のモニター環境の聴取位置での周波数特性を調べたい場合、マイクを聴取位置に設置し、スピーカーから測定信号を放出し、信号を収録することでインパルス応答が得られます。得られたインパルス応答には、室内が時不変だと仮定すれば「出力機器→室内音場→収録機器」の情報が全て含まれているということになります。

また、こういうおカタイ分析だけでなく、色々な場所で測定したインパルス応答を任意の音源に畳み込めばトリッキーなリバーブを作れたりします。

例えば、The Open Acoustic Impulse Response Library: OpenAIRには教会やトンネルの中などの実空間のインパルス応答が公開されてます。リバーブプラグインの代表格のAltiverbも実空間のインパルス応答を使用したコンボリューションリバーブです。

今回はPythonでインパルス応答測定を行うプログラムを作成しました。(実は私の指導教官が作ったMatlabスクリプトを移植しただけでほぼ丸写しなんですが、良いPython練習になりました。)

https://github.com/imaimamu/pinktspAnalyzer

pinktspGenerator.pyでPink-TSP信号を生成する

pinktspGenerator.pyは、インパルス応答の測定法の一つであるPink-TSP法(Pink-Time-Stretched-Pulse method)で使用するPink-TSP信号を生成し、wavで書き出すものです。測定に使う信号についてはかなり色々な考え方があって日々改善されているので、Pink-TSPが全てではないですが、割とポピュラーだと言われています。

前回の記事で作ったnextpow2と、素数周りの計算をしてくれるprimes.py、wavデータでの書き出しにutility.pyが必要になります。

"""
pinktspGenerator.py
Generate pinktsp signal.
This code was originaly written by Atsushi Marui in Matlab,
and ported to Python by Hidetaka Imamura.

Copyright 2009 MARUI Atsushi
Copyright 2015 Hidetaka Imamura
"""

import nextpow2
import primes
import cmath
import numpy as np
from scipy.io import wavfile
from utility import float2pcm

# User settings
dur = 5      # length of signal (seconds)
fs = 48000   # number of samples per second (Hz)
nbits = 16   # number of bits per sample (bit)
reps = 4     # number of repeated measurements (times)


N = 2 ** (nextpow2.nextpow2(dur * fs))
m = primes.primes_below(N / 4)
m = m[-1]

a = 2 * m * np.pi / ((N / 2) * np.log(N / 2))
j = cmath.sqrt(-1)
pi = (cmath.log(-1)).imag

H = np.array([1])
H = np.hstack([H, np.exp(j * a * np.arange(1, N / 2 + 1) * np.log(np.arange(1, N / 2 + 1))) / np.sqrt(np.arange(1, N / 2 + 1))])
H = np.hstack([H, np.conj(H[(N / 2 - 1):0:-1])])
h = (np.fft.ifft(H)).real
mvBefore = np.abs(h)
mv = min(mvBefore)
mi = np.where(mvBefore == mvBefore.min())
mi = mi[0]
h = np.hstack([h[mi:len(h)], h[0:mi]])
h = h[::-1]

Hinv = np.array([1])
Hinv = np.hstack([Hinv, np.exp(j * a * np.arange(1, N / 2 + 1) * np.log(np.arange(1, N / 2 + 1))) * np.sqrt(np.arange(1, N / 2 + 1))])
Hinv = np.hstack([Hinv, np.conj(Hinv[(N / 2 - 1):0:-1])])
hinv = (np.fft.ifft(Hinv)).real

hh = np.hstack((np.tile(h, (reps, 1)).flatten(), np.zeros(len(h))))
out = hh / max(np.abs(hh)) / np.sqrt(2)

wavfile.write('pinktsp.wav', fs, float2pcm(out, 'int16'))

 

こんな波形のWAVデータが出力されました。(図はAudacityで作成しています)

スクリーンショット 2015-07-10 22.46.19

 

※音量に注意して下さい。

 

できあがったPink-TSP信号を使って測定、分析する方法はまた次回の投稿で!

This Post Has 0 Comments

Leave A Reply