pythonでインパルス応答測定~直線状畳み込みver.~

今回はこちらの記事で紹介したインパルス応答測定の直線状畳み込みを用いた手法を紹介します。

円状畳み込みを用いた手法の実装は以下の記事からチェック可能です。

pythonでインパルス応答測定~円状畳み込みver.~

インパルス応答の概要も上記の記事で解説しているので参考にしてください。

記事の参考資料は引き続きこちらを参考にしています。

直線状畳み込みの原理を用いたインパルス応答測定のメリットは以下にあります。

• 用いる信号長とインパルス応答の長さに大きな制限がなく、インパルス応答長よりも短い信号で測定できる
• 上の利点により短めの信号で実際に音場の響きを聞きながら測定できる

じゃあ、前回の円状畳み込みの原理を用いた手法のメリットはあるのか？と考えましたが、

インパルス応答が短いときにはこちらでやる必要はないんじゃないかと思います。

では、早速実装に移りたいと思いますが

基本的には測定信号を並べて同期加算して逆フィルタという流れです。

円状畳み込みとの違いは

• 測定信号に間隔を空けて並べる（空ける間隔はインパルス応答長よりも長く設定する）
• 逆フィルタとの畳み込みに直線状畳み込みを行う

これだけです。前回の記事と比較しながら実装すると分かりやすいと思います。

今回は測定したいインパルス応答長をサンプリングレートの3000ポイント分として、信号長を2048としました

まずは以下をインポート

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as fft
import soundfile as sf
import math

Swept-Sine信号を作成していきます。

def Swept-Sine():
    gain = 70
    TSP = 2048#信号長
    m_exp = 1 
    m = TSP // (2 ** m_exp)  # (J=2m)実効長m=(3/4)N~(1/2)N程度
    a = ((m * np.pi) * (2.0 / TSP) ** 2)#α=4mπ/N^2
    ss_freqs = np.zeros(TSP, dtype=np.complex128)
    ss_freqs[:(TSP // 2) + 1] = np.exp(-1j * a * (np.arange((TSP // 2) + 1) ** 2))#exp(jαk^2)
    ss_freqs[(TSP // 2) + 1:] = np.conj(ss_freqs[1:(TSP // 2)][::-1])
    # ifft and real
    ss = np.real(fft.ifft(ss_freqs))#逆フーリエ変換で信号作成
    # roll
    ss = gain * np.roll(ss, TSP//2 - m)#信号をシフトさせる
    return(ss)

これでサンプリングレート\(f_s\)の1/2まで遷移する信号が作成できます。（保存時や再生時にサンプリングレートを指定する）
このssをwavで保存するか、呼び出すかして再生できるようにしておきます。
上記で作った信号ssを用いて、インパルス応答の測定のための準備をしていきます。

まずは測定信号作成。小文字が時間領域、大文字が周波数領域であることを意識して読んでください。

"""測定信号作成"""
ss = np.r_[ss, np.zeros(3500)]#3500ポイント間隔を空ける
ss_rep =  np.tile(ss, 3)#同期加算回数分つなげる

ここが前回との違いの一つ目です。

測定信号の間に間隔を空けましょう。今回は3000+500ポイント空けました。
次に逆フィルタを作成（DFTでの直線畳み込みのため０パディングをします）

inv_SS = np.conj(fft.fft(ss))#周波数領域で複素共役
inv_SS =  np.r_[inv_SS, np.zeros(3500)]

0パディングして、直線畳み込みを実現します。

あとはss_repを再生して、録音します。自由にやってもらっても構いませんが、信号長を厳密に切り出すには同時に録音と再生を行う必要があります。

オーディオインターフェースを用いて同時録音・再生を行う方法は以下の記事で解説しています。

ASIOなオーディオデバイスで同時録音・再生(python)

pythonにて音声を同時に再生・録音する方法をまとめました。sounddeviceというモジュールを使って実現しました.KOMPLETE AUDIO6をwindowsで使うときはASIOに対応させる必要があるのでそのあたりを説明しています。

得た出力をoutputとします。

加算平均を行って、フーリエ変換します。

加算平均で取り込む長さは、信号長＋間隔長（今回は2048+3500）

"""出力の加算平均"""
start = 0
output_ave = np.zeros(TSP)
for i in range(3):並べた信号の数
    output_ave += output[start:start+TSP+3500]
    start += TSP+3500
output_ave /= 3#加算したものを割る（加算平均）

得た加算平均のフーリエ変換と逆フィルタの線形畳み込みを行います。

"""伝達関数TF"""
TF = inv_SS*fft.fft(output_ave)
"""インパルス応答IR"""
IR = fft.ifft(TF)

以上でインパルス応答を測定することができます。
みなさんも試してみてください。