Python：マルチバンド処理用フィルタバンクのクラス設計・実装・動作検証

以前の記事では，木構造のような信号処理フローのマルチバンド処理を紹介しました．

この記事では、実装・利用しやすさに着目して信号処理フローを再設計・実装します。

設計

以前の記事では木構造の処理フローでしたが，今回は低域のカットオフ周波数から徐々に算出するような構造とします。

オールパスフィルタの伝達関数をまとめると以下のようになり，各帯域の信号はそれぞれの最後のオールパスフィルタで位相が同期されることで、回路全体が元信号の振幅成分を変化させないオールパスフィルタ回路となります。

なお、各オールパスフィルタのフィルタ係数はこちらの記事にあるようにローパス・ハイパスフィルタのフィルタ係数から簡単に算出できます．

www.wizard-notes.com

実装

フィルタバンククラスの設計

マルチバンド処理用フィルタバンクのクラスを作ることで，使い勝手のよい実装にします．

クラス定義・コンストラクタは以下の通りです．

今回は簡単化のため，奇数次バターワースのみ利用可能とします．

class AudioCrossoverFilters():
    def __init__(
            self,
            sr, 
            order=3,
            worN=10000,
        ):
        self.sr = sr
        self.order = order
        self.worN = worN

        if self.order % 2 == 0:
            sys.exit()

フィルタ係数・伝達関数算出部

フィルタ係数・伝達関数算出部では，フィルタバンクのフィルタ数に応じてローパス・ハイパス・オールパスフィルタ係数を算出します．

以下の記事より，バターワースフィルタで補完フィルタを作成する場合は，分母項（自己回帰フィルタ）のフィルタ係数aを共通化することができます．

www.wizard-notes.com

def compute_complementary_filters(self, fc):
    # フィルタ係数・伝達関数算出(バターワースフィルタ)
    b_lp, a = scipy.signal.butter(self.order, fc/(self.sr/2), btype='low')
    b_hp, _ = scipy.signal.butter(self.order, fc/(self.sr/2), btype='high')
    b_ap    = b_lp + b_hp
    b       = [b_lp, b_hp, b_ap]
   
    # 伝達関数
    h = [None] * 3
    w, h[0] = scipy.signal.freqz(b_lp, a, worN=self.worN)
    _, h[1] = scipy.signal.freqz(b_hp, a, worN=self.worN)
    _, h[2] = scipy.signal.freqz(b_ap, a, worN=self.worN)
    
    self.w = w
    
    return b, a, h

伝達関数算出

フィルタバンクの処理フローに従って伝達関数を計算します．

以下の関数では低域のフィルタから順番に算出しています．

def generate(self, fc):
    # カットオフ周波数 fc を numpy.array化
    if type(fc) is not list and type(fc).__module__ != np.__name__:
        fc = [fc]
    self.fc_list = np.array(fc)
    self.fc_list.sort()
    print(fc)
    
    self.n_filters = len(self.fc_list)
    self.n_bands   = len(self.fc_list) + 1
    
    # フィルタ係数・伝達関数の格納リスト
    self.b_lp_list = [None] * self.n_filters
    self.b_hp_list = [None] * self.n_filters 
    self.b_ap_list = [None] * self.n_filters
    self.a_list    = [None] * self.n_filters
    self.h_lp_list = [None] * self.n_filters
    self.h_hp_list = [None] * self.n_filters
    self.h_ap_list = [None] * self.n_filters

    for k, fc in enumerate(self.fc_list):
        b, a, h = self.compute_complementary_filters(fc)
        self.b_lp_list[k] = b[0]
        self.b_hp_list[k] = b[1]
        self.b_ap_list[k] = b[2]            
        self.a_list[k]    = a
        self.h_lp_list[k] = h[0]
        self.h_hp_list[k] = h[1]
        self.h_ap_list[k] = h[2] 

    # 回路全体の伝達関数を計算
    self.h_all = 0.0 
    for k in range(self.n_bands):
        h_tmp = 1.0
        if k < self.n_filters:
            h_tmp *= self.h_lp_list[k] # lowpass
            #print(f"{k}: add lowpass[{k}]")
            
        for m in range(0, k):
            h_tmp *= self.h_hp_list[m] # highpass
            #print(f"{k}: add highpass[{m}]")
            
        for m in range(k+1, self.n_filters):
            h_tmp *= self.h_ap_list[m] # allpass
            #print(f"{k}: add allpass[{m}] {self.n_filters}")
        
        self.h_all += h_tmp 

フィルタリング処理

入力信号から各帯域の信号を算出します．

こちらも先ほどと同様に低域から算出しています．

なお，信号再構成の検証における入出力差分計算時の遅延合わせが面倒なのでゼロ位相特性フィルタを使っていますが，通常のフィルタリング (scipy.signal.lfilter)でも構いません．

www.wizard-notes.com

def filt(self, x):
    tmp_x = x.copy()
    y = np.zeros((self.n_bands, x.shape[0]))
    filt = scipy.signal.filtfilt #scipy.signal.lfilter

    for k in range(self.n_filters):
        # Lowpass
        b, a = self.b_lp_list[k], self.a_list[k]
        y[k,:] = filt(b, a, tmp_x)
        # Allpass
        for m in range(k+1, self.n_filters):
            b, a = self.b_ap_list[m], self.a_list[m]
            y[k,:] = filt(b, a, y[k])
             
        b, a = self.b_hp_list[k], self.a_list[k]
        tmp_x = filt(b, a, tmp_x)
    y[self.n_bands-1, :] = tmp_x
        
    return y

検証結果

周波数応答のプロット

以下のようなスクリプトによって周波数応答をプロットすることができます．

周波数応答のプロットに使っているbode_plotはこちらの記事を参照してください．

www.wizard-notes.com

sr = 44100
order = 3
n_bands = 9
start = 2 # 最初のカットオフ周波数を10**startで指定
resolv = 0.25 # オクターブの分割指定。0.25なら4分割
fc = 10 ** (resolv * np.arange(0, n_bands-1) + start) 

acf = AudioCrossoverFilters(sr, order=order)  
acf.generate(fc)

# plot
bode_plot(
    acf.w, 
    acf.h_all, # 回路全体の周波数応答のみプロット
    #acf.h_lp_list + acf.h_hp_list + [acf.h_all], #全てのフィルタの周波数応答をプロット
    #labels=labels,
    sr=sr,
    fc=fc,
    xscale="log"
)

2バンド（カットオフ周波数1つ）

n_bands = 2の時の結果です

4バンド（カットオフ周波数3つ）

8バンド

信号へのフィルタリング

1秒のチャープ信号を使って目検で確かめると，低域から順にそれぞれの帯域の信号が分配され，また，再合成できることを確認しました．

import soundfile as sf
x, sr = sf.read("short_chirp.wav")
y = acf.filt(x)
# y.shape = (バンド数, 信号サンプル数)

y_sum = np.zeros(x.shape)
for k in range(y.shape[0]):
    y_sum += y[k]

以下は4バンドの時の結果です．

まとめ

マルチバンド処理用フィルタバンクのクラス設計・実装・動作検証をしました．

注意点としては，バターワースフィルタの次数は大きくても 3 か 5 次くらいがよいです。

それ以上大きくなると，量子化誤差の影響でフィルタが不安定になる傾向を確認しています．

なお，補完フィルタの性質として，H_lp_n+H_hp_n=H_ap_n という等式が成立するため，例えば　H_ap_n=H_ap_n-H_hp_n として片方が減算である並列接続を利用した回路も考えられます．

コード全体

import sys
import numpy as np
import scipy.signal

from bode_plot import bode_plot 

class AudioCrossoverFilters():
    def __init__(
            self,
            sr, 
            order=3,
            worN=10000,
        ):
        self.sr = sr
        self.order = order
        self.worN = worN

        if self.order % 2 == 0:
            sys.exit()
            
    def compute_complementary_filters(self, fc):
        # フィルタ係数・伝達関数算出(バターワースフィルタ)
        b_lp, a = scipy.signal.butter(self.order, fc/(self.sr/2), btype='low')
        b_hp, _ = scipy.signal.butter(self.order, fc/(self.sr/2), btype='high')
        b_ap    = b_lp + b_hp
        b       = [b_lp, b_hp, b_ap]
       
        # 伝達関数
        h = [None] * 3
        w, h[0] = scipy.signal.freqz(b_lp, a, worN=self.worN)
        _, h[1] = scipy.signal.freqz(b_hp, a, worN=self.worN)
        _, h[2] = scipy.signal.freqz(b_ap, a, worN=self.worN)
        
        self.w = w
        
        return b, a, h
        
    def generate(self, fc):
        # カットオフ周波数 fc を numpy.array化
        if type(fc) is not list and type(fc).__module__ != np.__name__:
            fc = [fc]
        self.fc_list = np.array(fc)
        self.fc_list.sort()
        print(fc)
        
        self.n_filters = len(self.fc_list)
        self.n_bands   = len(self.fc_list) + 1
        
        # フィルタ係数・伝達関数の格納リスト
        self.b_lp_list = [None] * self.n_filters
        self.b_hp_list = [None] * self.n_filters 
        self.b_ap_list = [None] * self.n_filters
        self.a_list    = [None] * self.n_filters
        self.h_lp_list = [None] * self.n_filters
        self.h_hp_list = [None] * self.n_filters
        self.h_ap_list = [None] * self.n_filters

        for k, fc in enumerate(self.fc_list):
            b, a, h = self.compute_complementary_filters(fc)
            self.b_lp_list[k] = b[0]
            self.b_hp_list[k] = b[1]
            self.b_ap_list[k] = b[2]            
            self.a_list[k]    = a
            self.h_lp_list[k] = h[0]
            self.h_hp_list[k] = h[1]
            self.h_ap_list[k] = h[2] 

        # 回路全体の伝達関数を計算
        self.h_all = 0.0 
        for k in range(self.n_bands):
            h_tmp = 1.0
            if k < self.n_filters:
                h_tmp *= self.h_lp_list[k] # lowpass
                #print(f"{k}: add lowpass[{k}]")
                
            for m in range(0, k):
                h_tmp *= self.h_hp_list[m] # highpass
                #print(f"{k}: add highpass[{m}]")
                
            for m in range(k+1, self.n_filters):
                h_tmp *= self.h_ap_list[m] # allpass
                #print(f"{k}: add allpass[{m}] {self.n_filters}")
            
            self.h_all += h_tmp

    def filt(self, x):
        tmp_x = x.copy()
        y = np.zeros((self.n_bands, x.shape[0]))
        filt = scipy.signal.filtfilt #scipy.signal.lfilter

        for k in range(self.n_filters):
            # Lowpass
            b, a = self.b_lp_list[k], self.a_list[k]
            y[k,:] = filt(b, a, tmp_x)
            # Allpass
            for m in range(k+1, self.n_filters):
                b, a = self.b_ap_list[m], self.a_list[m]
                y[k,:] = filt(b, a, y[k])
                 
            b, a = self.b_hp_list[k], self.a_list[k]
            tmp_x = filt(b, a, tmp_x)
        y[self.n_bands-1, :] = tmp_x
            
        return y


sr = 44100
order = 7
n_bands = 8
start = 2 # 最初のカットオフ周波数を10**startで指定
resolv = 0.25 # オクターブの分割指定。0.25なら4分割
fc = 10 ** (resolv * np.arange(0, n_bands-1) + start) 

acf = AudioCrossoverFilters(sr, order=order)  
acf.generate(fc)

# plot
bode_plot(
    acf.w, 
    #acf.h_all, # 回路全体の周波数応答のみプロット
    acf.h_lp_list + acf.h_hp_list + [acf.h_all], #全てのフィルタの周波数応答をプロット
    #labels=labels,
    sr=sr,
    fc=fc,
    xscale="log"
)

import soundfile as sf
x, sr = sf.read("shortchirp.wav")
y = acf.filt(x)
# y.shape = (バンド数, 信号サンプル数)

y_sum = np.zeros(x.shape)
for k in range(y.shape[0]):
    y_sum += y[k]