【python】頭部伝達関数を重畳保留法を用いて畳み込んでみた

今回は重畳保留法(overlap_save)を使って頭部伝達関数を畳み込んでみました。

・頭部伝達関数

前回と同じ

・重畳保留法

以下のサイトを参考に書きました。
原理はわかっていません…
元気があったら勉強してみようと思います

www.softist.com

上記のサイトを参考にパラメータを決めていくと

N(fftブロック長) = 1024
L(フィルタ長) = 512
M(音源の切り出し区間) = 1024

MとNが同じなのでプログラムではMの部分にNを使いました。

・ソースコード

今回も平面のみでやりました

プログラムの流れは以下のとおりです

・音源からM点抜き出す
・L点のhrtfの左側を0で埋めてN点の配列を作る
・この2つをfftして掛けあわせ、ifftする
・ifftして得られたデータの最初からL点を抜き出してストリームに流し込む

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

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(soundData, hrtf, N, L):
    tmpFilter = numpy.zeros(N)
    tmpFilter[L:] += hrtf

    spectrum = numpy.fft.fft(soundData, n=N)
    hrtfFft = numpy.fft.fft(tmpFilter, n=N)
    add = spectrum * hrtfFft
    result = numpy.fft.ifft(add, n=N).real
    return result[:L]

def play_elev0(soundData, hrtfL, hrtfR, L, N, streamObj, position):
    index = 0
    while(soundData[index:].size > L):
        resultData = numpy.empty((0, 2), dtype=numpy.int16)

        convL = convolution(soundData[index:index + N, 0], hrtfL[position], N, L)
        convR = convolution(soundData[index:index + N, 1], hrtfR[position], N, L)

        for i in range(convL.size):
            resultData = numpy.append(resultData, numpy.array([[int(convL[i]), int(convR[i])]], dtype=numpy.int16), axis=0)

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

    streamObj.close()

#wavファイルのサウンドデータをロード
soundDataPath = "ongaku.wav"
rate, soundData = wav.read(soundDataPath)

#HRTFのロード
hrtfL, hrtfR = load_elev0hrtf()

#PyAudioの初期設定
p = pyaudio.PyAudio()
stream = p.open(format = 8,
                channels = 2,
                rate = rate,
                output = True)
#パラメータ設定
L = 512
N = 1024
#音の方向(左回り)
position = 18

play_elev0(soundData, hrtfL, hrtfR, L, N, stream, position)
p.terminate()