Pythonでスピーカの周波数特性を分析

なぜスピーカの周波数分析を行うのか

ダンスミュージックを聴くときは低音重視、クラシックを聴くときは高音重視、音楽制作に使うときはフラットなど、スピーカは用途によって求められる周波数の特性が変わってきます。欲しい音がきちんと出るスピーカを設計したり、選んだりする際にスピーカの周波数分析を行い、定量的に評価することが重要になります。

分析の流れ

PythonでTSP信号（測定信号）を作る

→スピーカに入力して出力をマイクで録音

→Pythonで分析、グラフ作成

使用機材

• 被測定スピーカ：Genelec 8050A
• マイク：NTI M2230
• オーディオインターフェイス：RME MADIface XT
• 測定環境：無響室

TSP信号とは

システムの周波数特性は、理想的にはインパルスをシステムに入力し、その出力であるインパルス応答をフーリエ変換することで得られます。
しかしインパルスを放出するスピーカの出力が十分でないなどの条件があり一般的にはこれを実際に行うのは難しいです。
TSP(Time Stretched Pulse)信号は直接的にインパルスを使わずにインパルス応答を算出するときに用いる信号であり、インパルスを時間的に引き伸ばしたような信号で、周波数領域において以下の式で与えられます。

S(k) = \left\{
\begin{array}{ll}
\exp \Bigl(\frac{-j4m\pi k^2}{N^2} \Bigr) & (0 \leq k \leq \frac{N}{2}) \\
S^*(N - k) & (\frac{N}{2} \lt k \lt N)
\end{array}
\right.

ここで　$S^*(k)$は$S(k)$の複素共役、 $m = N/4, N = 2^l(lは整数)$, であり、NはFFTの点数です。
これを逆フーリエ変換することによって被測定スピーカに入力するためのTSP信号の時間波形$s[n]$を得ることができます。
このTSP信号を無響室内で被測定スピーカで再生しマイクで録音すると出力信号$y[n]$が得られ、$y[n]$のフーリエ変換を$Y(k)$、求める周波数特性を$H(k)$とすると

\begin{align}\\
Y(k) = H(k)S(k)\\
\therefore H(k) = \frac{Y(k)}{S(k)}\\
 \Bigl(=Y(k)S^{-1}(k) \Bigr)

\end{align}

となりTSP信号の逆特性をもつ信号のフーリエ変換と、出力信号のフーリエ変換の積を計算することによってスピーカーの周波数特性を求めることができます。
ここで$S^{-1}(k)$は

S^{-1}(k) = \left\{
\begin{array}{ll}
\exp \Bigl(\frac{j4m\pi k^2}{N^2} \Bigr) & (0 \leq k \leq \frac{N}{2}) \\
S^{-1*}(N - k) & (\frac{N}{2} \lt k \lt N)
\end{array}
\right.

で求められます。
またTSP信号のバリエーションとしてPink TSPというものも存在し、こちらは高音域でのS/N比が悪くなる代わりに広帯域のバランスが取れるという特徴を持っています。
Pink TSPの周波数特性は以下の式で与えられます。

S(k) = \left\{
\begin{array}{ll}
1 & (k=1)\\
\frac{\exp(jak\log{k})}{\sqrt{k}} & (0 \lt k \leq \frac{N}{2}) \\
S^*(N - k) & (\frac{N}{2} \lt k \lt N)
\end{array}
\right.\\
S^{-1}(k) = \left\{
\begin{array}{ll}
1 & (k=1)\\
\sqrt{k}\exp(jak\log{k}) & (0 \lt k \leq \frac{N}{2}) \\
S^{-1*}(N - k) & (\frac{N}{2} \lt k \lt N)
\end{array}
\right.

今回はTSP信号をPythonで作っていきます。使うライブラリはNumpy、SoundFileです。

tsp_design.py

import numpy as np

def normal_tsp(n, gain=50, repeat=1):
    N = 2**n
    m = N//4

    A = 50
    L = N//2 - m
    k = np.arange(0, N)

    tsp_freq = np.zeros(N, dtype=np.complex128)
    tsp_exp = np.exp(-1j*4*m*np.pi*(k/N)**2)

    tsp_freq[0:N//2+1] = tsp_exp[0:N//2+1]
    tsp_freq[N//2+1: N+1] = np.conj(tsp_exp[1 : N//2][::-1])

    tsp_inv_freq = 1 / tsp_freq

    tsp = np.real(np.fft.ifft(tsp_freq))
    tsp = gain * np.roll(tsp, L)

    tsp_repeat = np.r_[np.tile(tsp, repeat), np.zeros(N)]

    tsp_inv = np.real(np.fft.ifft(tsp_inv_freq))
    tsp_inv =  gain * np.roll(tsp_inv, -L)

    tsp_inv_repeat = np.r_[np.tile(tsp_inv, repeat), np.zeros(N)]

    return tsp_repeat, tsp_inv

def pink_tsp(n, gain=50, repeat=1):

    N = 2**n
    m = N//4

    L = N//2 - m
    k = np.arange(1, N)

    a = 4 * m * np.pi / (N * np.log(N/2))

    tsp_freq = np.zeros(N, dtype=np.complex128)
    tsp_exp = np.exp(1.j * a * k * np.log(k)) / np.sqrt(k)

    tsp_freq[0] = 1
    tsp_freq[1:N//2+1] = tsp_exp[1:N//2+1]
    tsp_freq[N//2+1: N+1] = np.conj(tsp_exp[1 : N//2][::-1])

    tsp_inv_freq = 1 / tsp_freq

    tsp = gain * np.real(np.fft.ifft(tsp_freq))[::-1]
    tsp =  gain * np.roll(tsp, L)
    tsp_repeat = np.r_[np.tile(tsp, repeat), np.zeros(N)]

    tsp_inv = np.real(np.fft.ifft(tsp_inv_freq))[::-1]
    tsp_inv =  gain * np.roll(tsp_inv, L)


    return tsp_repeat, tsp_inv

このtsp_design.pyとSoundFileをインポートして、関数を用いて生成された配列をwavに変換します。
今回は出力信号のS/N比を高くするために加算平均を取りたいので、同じ信号を５回繰り返させます。

>>> import tsp_design
>>> import soundfile as sf
>>> normal_tsp, normal_tsp_inv = tsp_design.normal_tsp(18, repeat=5)
>>> pink_tsp, pink_tsp_inv = tsp_design.pink_tsp(18, repeat=5)
>>> sf.write('normal_tsp.wav', normal_tsp, 48000)
>>> sf.write('pink_tsp.wav', pink_tsp, 48000)

これで出来上がりです。

測定

測定のブロックダイアグラムは以下のようになります。

PC上のDAW（今回はLogic Pro X）から生成したTSP信号をスピーカに入力し、出力信号をマイクで録音していきます。

分析

今回はS/N比を向上させるために同期加算を行います。具体的にはスピーカからの出力信号を入力したTSP信号の長さだけ切り分けそれらの加算平均を取ることによってノイズ成分を打ち消します。


この出力波形のフーリエ変換とTSP信号の逆特性の積を計算するとスピーカの周波数特性が得られます。

ちなみに、この周波数特性を逆フーリエ変換すると、システムのインパルス応答が得られます。

以上の一連の処理は以下のコードで行なえます。

analysis.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

from tsp_design import normal_tsp, pink_tsp
import soundfile as sf

def sychronous_addition(filename, repeat, N):
    '''
    filename : Name of wav file (str)
    repeat : Number of repetition (int)
    N : Length of input signal (int)
    '''
    data, fs = sf.read(filename)

    # add zeros if length is too short
    if len(data) < repeat * N:
        data = np.r_[data, np.zeros(repeat * N - len(data))]

    mean = np.zeros(N)
    for i in range(repeat):
        mean += data[i * N : (i + 1) * N]
    mean = mean / repeat

    return mean

n = 18
N = 2 ** n
fs = 48000

normal_mean = sychronous_addition('normal_tsp_response.wav', 5, N)
pink_mean = sychronous_addition('pink_tsp_response.wav', 5, N)
normal_tsp, normal_tsp_inv = normal_tsp(n)
normal_tsp_inv_freq = np.fft.fft(normal_tsp_inv)

pink_tsp, pink_tsp_inv = pink_tsp(n)
pink_tsp_inv_freq = np.fft.fft(pink_tsp_inv)

H_normal = np.fft.fft(normal_mean) * normal_tsp_inv_freq
H_pink = np.fft.fft(pink_mean) * pink_tsp_inv_freq

h_normal = np.fft.ifft(H_normal)
h_pink = np.fft.ifft(H_pink)

f = np.linspace(0, fs/2, N//2)

plt.plot(f[1:], 20*np.log10(np.abs(H_normal[1:N//2])/np.max(np.abs(H_normal))), color="red", label="Nornal TSP")
plt.plot(f[1:], 20*np.log10(np.abs(H_pink[1:N//2])/np.max(np.abs(H_pink))), color="blue", label="Pink TSP")
plt.xscale('log')
plt.xlim(20, 20000)
plt.title('Genelec 8050Aの周波数特性')
plt.xlabel('周波数[Hz]')
plt.ylabel('レベル[dB]')
plt.legend()
plt.show()

plt.plot(h_normal)
plt.title("Inpulse Response(with normal TSP)")
plt.xlim(0, 10000)
plt.show()

plt.plot(h_pink)
plt.title("Inpulse Response(with Pink TSP)")
plt.xlim(0, 10000)
plt.show()

まとめ・考察

今回は2種類のTSP信号を用いてスピーカの周波数特性を分析してみました。
結果のグラフから見ると、Pink TSPで測定したほうが100Hz以下の低域のスペクトルがきれいに出ています。
また次回からいろいろな測定手法も検討していきたいと思います。