cBlog

Tips for you.

【測定】Logic Pro付属のリニアフェイズEQは厳密に線形位相特性

スポンサーリンク

リニアフェイズEQで位相が乱れるはずがないと思っていたのですが、もしかして近似線形位相ってこともあるのか?と思い、公称を鵜呑みにしている自分に気がついたので測定してみました。

対象はLogic Pro付属の「Linear Phase EQ」。これしかインストールしていないんだもの。結果はちゃんと線形位相特性でした。というわけで、位相乱れ・波形歪みはありません。

 

測定方法

量子化以外誤差の混入し得ないデジタル領域の測定なので定義通りインパルス信号を入力してやってもよい(はずな)のですが、定石通りサインスイープ(TSP、Swept Sine)信号を用いて測定しました。サインスイープ信号は以下のような振幅が1、周波数が時間に比例する(位相が周波数の2乗に比例する)信号です。

\begin{align}
S(k) =
\begin{cases}
e^{-j2\pi J(k/N)^2}, &k = 0, 1,\ldots, \frac{N}{2} \\
S(N - k)^*, &k = \frac{N}{2} + 1, \frac{N}{2} + 2,\ldots, N - 1
\end{cases}
\end{align}

詳しくは以下の金田先生の講義資料を見てください。というか、これを見れば今回の測定ぜんぶできます。

http://www.asp.c.dendai.ac.jp/ASP/IRseminor2016.pdf

Logic Proでは、ポインタを置くとプラグインの遅延量が確認できます。Linear Phase EQでは「2560 個のサンプル」のように表示されるので、ここから群遅延が2560、フィルタ次数5120、インパルス応答長5121と推察できます。そこで、信号長N=8192、実効長J=4096としてサインスイープ信号を生成しました。また、立ち上がりを滑らかにするためにサインスイープ信号は2048サンプルでロールしてあります。詳しくは金田先生の(ry

S^(-1)(k)は安定で簡単に求まるので、出力信号をDFT(実際はFFTだがFFTはアルゴリズム名)したものに掛けてやれば周波数特性H(k)を求めることができます。

補足:一般にインパルス応答長は無限か読めないので、出力信号(長さ:入力信号長+インパルス応答長-1)はDFT長をはみ出す可能性がある。DFTの周期性により、畳み込み長がDFT長を超える場合、DFT領域の積は循環畳み込みとなる。それに合わせ、線形畳み込みである被測定系で循環畳み込みを行うために、2周期再生し、出力の2周期目を使う。

実際にLogic Proでリージョンをオーディオファイルとして書き出し…たのですが、遅れない… 環境設定 > オーディオ > 一般 > プラグインのレイテンシの「プリロールを再生」をオフにしたらちゃんと2560サンプル遅れました。プラグインの遅延量早めて再生する機能のようです。念のため、同階層にある「補正」もオフにして計測しています。

 

測定

手始めに500 Hzだけを3dBブースト。

Linear Phase EQ設定 500 Hz +3dB

位相特性はこちら(単位はrad)。

Linear Phase EQ位相特性 500 Hz +3dB

直線ですね。次に、群遅延で遅延量を見てみます。

群遅延とは、周波数に対してその周波数成分がどれだけ遅れるか表した量。角周波数での位相の微分の負。割り算だとその周波数までの平均遅延になるので、微分。位相だと振動がどの段階かを表すのですが、わかりづらい。高周波ほど、周期が長い低周波に合わせてめっちゃ位相を遅らせなきゃ波形が揃わないからです。群遅延だと単位はサンプル(詳細は後述)になるので便利です。まさに、群遅延は正弦波群の遅延を表しています。

こちらが群遅延特性。

Linear Phase EQ群遅延特性 500 Hz +3dB

群遅延は2560(縦軸の見方に注意)。Logic Proの表示通りでした。偏差も1サンプルと比べてずっと小さいので位相乱れはないと言っていいでしょう。

念のため、読み込み > 書き出ししただけの群遅延も求めてみました。

FXなし 群遅延特性

小さいですが同様に偏差はあるので、EQ時の群遅延の乱れも大方、インパルス応答算出時の演算誤差ではないでしょうか。f=0付近で暴れるのはよくわかりませんが、今回使った群遅延の数値計算の傾向のようです(起きないこともあるしf=fs/2で暴れることもある)。

余談ですが、サインスイープ信号は[-(2^23), 2^23-1]のフルスケールで生成したのですが、Logic Proに読み込んで書き出す(一度floatになる)と[-(2^23-1), 2^23-1]となるようです(ノーマライズをオフにしても)。2^(-23)だけダイナミックレンジ損すると思うのだが、それでいいのかLogicよ…

インパルス応答も出力してみました。

Linear Phase EQインパルス応答 500 Hz +3dB

[追記 2018/11/10] n=2560より前で値を持つことが、リニアフェイズEQのプリリンギングの原因となっています。

左右対称ですね。対称性を見るために、左右のRoot Mean Square Errorを計算してみました。結果は1.33*10^(-7)。24 bitの分解能2^(-23)と同オーダーですので、完全に対称と言えるはず。偶次偶対称の線形位相FIRフィルタですね。線形位相FIRフィルタの条件については以下の記事を。f=0、f=fs/2での特性の制約が異なります。

jp.mathworks.com

次は、プリセットの「06 Mastering > Final Mix - Pop」の群遅延。

Linear Phase EQ設定 Final Mix - Pop

Linear Phase EQ群遅延特性 Final Mix - Pop

EQ設定を変えても同様ですね。

比較のため、非リニアフェイズEQ「Channel EQ」でも以下の設定で群遅延を算出。

Channel EQ設定 500 Hz +3dB

Channel EQ群遅延特性 500 Hz +3dB

EQポイント前後で7サンプルくらい位相乱れが発生しています。線形位相じゃないとこうなるんです。

 

コラム:群遅延の単位

デジタル領域では、群遅延は無次元量であり、単位は無い(時間がサンプリング周期で正規化されていることを思い出してほしい)。あえて言えば、単位はサンプルとなる。[rad]を[rad/sample]で微分していると考えるとわかりやすいかもしれない。アナログ領域では、[rad]を[rad/s]で微分しているので単位は秒である。

 

ソースコード

(せっかくなので)

sine_sweep.py

サインスイープ信号を生成するプログラム。

時間領域でロールする代わりに、以下の性質を使って周波数領域でロールしています。

\begin{align}
f(n - n_0) \Leftrightarrow F(k)e^{-j2\pi n_0k/N}
\end{align}

import numpy as np
import matplotlib.pyplot as plt
import soundfile as sf


fs = 44100
N = 8192
J = N//2
roll = N//4

n = np.arange(N)
Y = np.empty(N, dtype=complex)
arg = -1j*2*np.pi*(J*(n[:N//2 + 1]/N)**2 + roll*n[:N//2 + 1]/N)
Y[:N//2 + 1] = np.exp(arg)
Y[N//2 + 1:] = np.conj(Y[N//2 - 1:0:-1])

y = np.real(np.fft.ifft(Y))
y /= np.max(np.abs(y))
y = np.tile(y, 2)

plt.plot(y)
plt.show()

sf.write('sine_sweep.wav', y, fs, 'PCM_24')

 

calc_gdelay.py

sine_sweep.pyで生成した信号の応答から群遅延特性をプロットするプログラム。

ここでは、群遅延D(k)の算出に以下のような位相のアンラッピングが不要なアルゴリズムを用いています。

\begin{align}
D(k) = \Re\left\{\frac{\mathrm{DFT}(k\cdot h(k))}{\mathrm{DFT}(h(k))}\right\}
\end{align}

導出は後で書く(かもしれない)。ざっくり言うと、H(k)の対数とると虚部が偏角になることを利用してD(k)を直接H(k)(とその導関数)で表して、インパルス応答をFIRフィルタ係数と考えたとき、微分すると肩の指数が係数として下りてくることを利用しています。

import sys
import numpy as np
import matplotlib.pyplot as plt
import soundfile as sf


if len(sys.argv) != 3:
    print('Usage: python {} in(sound) out(imresp)'.format(sys.argv[0]))
    sys.exit(1)

fni = sys.argv[1]
fno = sys.argv[2]

N = 8192
J = N//2
roll = N//4

n = np.arange(N)
X_inv = np.empty(N, dtype=complex)
arg = 1j*2*np.pi*(J*(n[:N//2 + 1]/N)**2 + roll*n[:N//2 + 1]/N)
X_inv[:N//2 + 1] = np.exp(arg)
X_inv[N//2 + 1:] = np.conj(X_inv[N//2 - 1:0:-1])

y, fs = sf.read(fni, start=N, frames=N)
Y = np.fft.fft(y)
H = Y*X_inv
phase = np.angle(H)
phase = np.unwrap(phase)

h = np.real(np.fft.ifft(H))
with open(fno, 'w') as flo:
    for i in range(N):
        print(h[i], file=flo)

gdelay = np.real(np.fft.fft(n*h)/H)

f = fs*n/N
#plt.plot(f[:N//2], phase[:N//2])
plt.plot(f[:N//2], gdelay[:N//2])
#plt.ylim(2559.96, 2560.04)
#plt.ylim(-0.05, 0.04)
plt.show()

 

次回予告

まず、今回の内容が前回の次回予告に適っていなかったことをおわびします。

次回は、標本化と量子化の影響についてお送りします。

アップしました

yaritakunai.hatenablog.com