はじめに

信号処理の分野では「最小位相フィルタ」というキーワードがでてきます。

その名の通り、フィルタのインパルス応答が時間的にできるだけ早く立ち上がり、最小限の遅延を持つように設計されたフィルタです。

Pythonのscipyライブラリには、便利な関数 scipy.signal.minimum_phase があり、これを使用することで簡単に最小位相フィルタを算出することができます*1。本記事では、この関数の基本的な使い方と内部で行われている計算の流れについて解説します。

どんな時に必要か

最小位相フィルタは、信号の時間的な遅延を最小限に抑える必要がある場合に使用されます。

例えば、特定の周波数特性を持つフィルタを設計し、それを最小位相フィルタに変換するケースを考えます。以下に、Pythonコードを使った算出・プロット例を示します。

クリックで展開

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal


# FIRフィルタの設計（窓関数法）
numtaps = 64  # フィルタのタップ数
cutoff = 0.3  # 正規化周波数でのカットオフ（0から0.5の範囲）
fir_filter = signal.firwin(numtaps, cutoff)

# 最小位相フィルタに変換
min_phase_filter = signal.minimum_phase(fir_filter, method='homomorphic')

# 周波数特性の計算
w, h_fir = signal.freqz(fir_filter, worN=1024)
w, h_min_phase = signal.freqz(min_phase_filter, worN=1024)

# 群遅延の計算
_, gd_fir = signal.group_delay((fir_filter, 1), w)
_, gd_min_phase = signal.group_delay((min_phase_filter, 1), w)
# プロット
plt.figure(figsize=(12, 8))

# インパルス応答
plt.subplot(3, 1, 1)
plt.plot(fir_filter, label='Original FIR Filter')
plt.plot(min_phase_filter, label='Minimum Phase FIR Filter', linestyle='--')
plt.title('Impulse Response')
plt.xlabel('Samples')
plt.ylabel('Amplitude')
plt.legend()
plt.grid()

# 周波数振幅特性
plt.subplot(3, 1, 2)
plt.plot(w / np.pi, 20 * np.log10(np.abs(h_fir)), label='Original FIR Filter')
plt.plot(w / np.pi, 20 * np.log10(np.abs(h_min_phase)), label='Minimum Phase FIR Filter', linestyle='--')
plt.title('Frequency Amplitude Response')
plt.xlabel('Normalized Frequency [×π rad/sample]')
plt.ylabel('Amplitude [dB]')
plt.legend()
plt.grid()

# 群遅延
plt.subplot(3, 1, 3)
plt.plot(w / np.pi, gd_fir, label='Original FIR Filter')
plt.plot(w / np.pi, gd_min_phase, label='Minimum Phase FIR Filter', linestyle='--')
plt.title('Group Delay')
plt.xlabel('Normalized Frequency [×π rad/sample]')
plt.ylabel('Group Delay [samples]')
plt.ylim([0,50])
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()

上記のコードでは、FIRフィルタを設計し、そのフィルタを最小位相に変換しました。プロットを見ると、最小位相フィルタのインパルス応答が早い時間に集中していることがわかります。また、振幅特性はほぼ同じで、群遅延は0に近づいていることが分かります。

つまり、最小位相フィルタ化することで、以下のメリットがあります。

フィルタによるシステムの遅延を削減
フィルタのタップ数削減
- 半分程度に減らすことができる

最小位相フィルタの算出方法

scipy.signal.minimum_phaseの紹介

scipy.signal.minimum_phase は、FIRフィルタのインパルス応答から最小位相のインパルス応答を生成するための関数です。

scipy.signal.minimum_phase(h, method='homomorphic')

h: フィルタのインパルス応答（通常はFIRフィルタの係数）。
method: 最小位相フィルタの計算方法。今回はデフォルトのモードである 'homomorphic' のみを使用します。

利用例

以下に、scipy.signal.minimum_phase を使用して最小位相フィルタを生成するPythonコードの例を示します。

import numpy as np
from scipy import signal

# FIRフィルタの係数を設定
numtaps = 64  # タップ数
cutoff = 0.3  # カットオフ周波数
fir_coeff = signal.firwin(numtaps, cutoff)

# 最小位相フィルタの生成
min_phase_fir = signal.minimum_phase(fir_coeff, method='homomorphic')

# 結果の表示
print("Original FIR Coefficients:\n", fir_coeff)
print("Minimum Phase FIR Coefficients:\n", min_phase_fir)

このコードは、まず通常のFIRフィルタを設計し、次に scipy.signal.minimum_phase 関数を使用して最小位相フィルタを生成します。method='homomorphic' を指定することで、ホモモルフィック処理（ヒルベルト変換による最小位相スペクトル推定）が適用され、最小位相の応答が得られます。

算出の流れ

ヒルベルト変換とホモモルフィック処理

scipy.signal.minimum_phase 関数での最小位相フィルタの算出には、ホモモルフィック処理（準同型フィルタリング）が使われます。これはヒルベルト変換による最小位相フィルタの推定となっています。

以下はそのコードと処理の流れです：

scipy/scipy/signal/_fir_filter_design.py at 92d2a8592782ee19a1161d0bf3fc2241ba78bb63 · scipy/scipy · GitHub

    else:  # method == 'homomorphic'
        # zero-pad; calculate the DFT
        h_temp = np.abs(fft(h, n_fft))
        # take 0.25*log(|H|**2) = 0.5*log(|H|)
        h_temp += 1e-7 * h_temp[h_temp > 0].min()  # don't let log blow up
        np.log(h_temp, out=h_temp)
        if half:  # halving of magnitude spectrum optional
            h_temp *= 0.5
        # IDFT
        h_temp = ifft(h_temp).real
        # multiply pointwise by the homomorphic filter
        # lmin[n] = 2u[n] - d[n]
        # i.e., double the positive frequencies and zero out the negative ones;
        # Oppenheim+Shafer 3rd ed p991 eq13.42b and p1004 fig13.7
        win = np.zeros(n_fft)
        win[0] = 1
        stop = n_fft // 2
        win[1:stop] = 2
        if n_fft % 2:
            win[stop] = 1
        h_temp *= win
        h_temp = ifft(np.exp(fft(h_temp)))
        h_minimum = h_temp.real

対数振幅スペクトル log|H| を計算
1を逆フーリエ変換
2の負の周波数を0埋め、正の周波数を2倍
3をフーリエ変換する
- 実部に元の振幅スペクトル，虚部に最小位相となる位相スペクトルが現れる
- log|H| + jφ(w) = log(|H| exp(jφ(w)))
4の指数関数をとった後、逆フーリエ変換して時間領域信号を得る
- exp(log(x)) = x
- 指数関数をとることで|H|exp(jφ(w))となるので、最小位相である信号が自然と得られる。