【python】頭部伝達関数を重畳加算法を使って畳み込んでみた

頭部伝達関数(hrtf)を重畳加算法(overlap_add)を使って音源に畳み込んでみました。

・頭部伝達関数

音源位置から両耳の入り口までの伝達関数

今回は以下のサイトで公開されているhrtfを使用させていただきました

・overlap_add

長い信号とfirフィルタの離散畳み込みを分割して効率的に処理する手法らしい

重畳加算法 - Wikipedia

ウィキを見てもらえばわかると思う

N = L + M - 1 = 2^p (N:fftブロック長、L:切り出し区間、M:フィルタ長)

と書かれている

使用するhrtfは512点なのでとりあえず

N = 1024

L = 513

M = 512

ウィキの図を見てもらうとオーバーラップ部分がM-1となっているので

overLap = 511

とした

・音源

フリーのwavファイル音源を使用しました

・ソースコード

とりあえず今回は水平面のhrtfのみ使用しました

import numpy
import pyaudio
import scipy.io.wavfile as scw

def load_elev0hrtf():
    elev0Hrtf_L = {}
    elev0Hrtf_R = {}

    for i in range(72):
        str_i = str(i * 5)

        if len(str_i) < 2:
            str_i = "00" + str_i
        elif len(str_i) < 3:
            str_i = "0" + str_i

        fileName = "L0e" + str_i + "a.dat"
        filePath = "../hrtfs/elev0/" + fileName
        test = open(filePath, "r").read().split("\n")

        data = []

        for item in test:
            if item != '':
                data.append(float(item))

        elev0Hrtf_L[i] = data

    for i in range(72):
        str_i = str(i * 5)

        if len(str_i) < 2:
            str_i = "00" + str_i
        elif len(str_i) < 3:
            str_i = "0" + str_i

        fileName = "R0e" + str_i + "a.dat"
        filePath = "../hrtfs/elev0/" + fileName
        test = open(filePath, "r").read().split("\n")

        data = []

        for item in test:
            if item != '':
                data.append(float(item))

        elev0Hrtf_R[i] = data

    return elev0Hrtf_L, elev0Hrtf_R

def convolution(data, hrtf, N, L):
    spectrum = numpy.fft.fft(data, n = N)
    hrtf_fft = numpy.fft.fft(hrtf, n = N)
    add = spectrum * hrtf_fft
    result = numpy.fft.ifft(add, n = N)
    return_data = result.real
    return return_data[:L], return_data[L:]

def play_elev0(sound_data, N, L, hrtfL, hrtfR, position, overLap, streamObj):
    index = 0
    overLap_L = numpy.zeros(overLap)
    overLap_R = numpy.zeros(overLap)

    while(sound_data[index:].size > L):
        result_data = numpy.empty((0, 2), dtype=numpy.int16)

        tmp_conv_L, add_L = convolution(sound_data[index:index + L, 0], hrtfL[position], N, L)
        tmp_conv_R, add_R = convolution(sound_data[index:index + L, 1], hrtfR[position], N, L)

        tmp_conv_L[:overLap] += overLap_L
        tmp_conv_R[:overLap] += overLap_R

        overLap_L = add_L
        overLap_R = add_R

        for i in range(tmp_conv_L.size):
            result_data = numpy.append(result_data, numpy.array([[int(tmp_conv_L[i]), int(tmp_conv_R[i])]], dtype=numpy.int16), axis=0)

        streamObj.write(bytes(result_data))
        index += L

    streamObj.close()


#初期設定など
soundDataPath = "./ongaku.wav"
rate, soundData = scw.read(soundDataPath)

p = pyaudio.PyAudio()
stream = p.open(format = 8,
                channels = 2,
                rate = rate,
                output = True)

hrtf_L, hrtf_R = load_elev0hrtf()
N = 1024
L = 513
overLap = 511
#音の再生する方向(0が正面で半時計回りで71まで)
position = 18

play_elev0(soundData, N, L, hrtf_L, hrtf_R, position, overLap, stream)
p.terminate()

hrtfをロードして畳み込んで再生といった流れです

間違いなどあればコメントください

zack_ff’s blog

ff14廃人がなにか書く

【python】頭部伝達関数を重畳加算法を使って畳み込んでみた