はじめに

マイクで収録した音にFIRフィルタで周波数補正をかけたくなりました。FIRフィルタというとscipy.signal.firwinなどがありますが、今回は単純なローパスフィルタやハイパスフィルタではなく、複雑な周波数特性の実現を目指します。私はFIRフィルタの理論を深く理解しているわけではないので、間違っている箇所があると思います。そのときはコメントで教えていただけるとありがたいです。

参考

東北大学講義資料の以下を参考にさせていただきました。
やる夫で学ぶディジタル信号処理

手順

ファイル読み込みなどのコードは省略しました。

1. マイク付属の周波数特性表の読み込み

これはとあるマイクの周波数特性と自由音場補正データです。

補正すべき値は以下のようになります。

2. 実現したい周波数特性をつくる

ここで、補正値は対数から線形になおしておきます。
このマイクは140 kHzまでが計測可能範囲なので、それより高い周波数ではゲインを0にします。

N = 128 #タップ数 + 1
fs = 500.0e3 #補正をかけたいデータのサンプリングレート
freq_array = np.linspace(0.0, fs, num=N)

gain_correction = np.zeros(N//2)
lim_freq = freq_array[freq_array <= 140.0e3]
# f_pro_funcとf_pro2_funcはスプライン補間テーブル
gain_correction[0:len(lim_freq)] = 10**(-(f_pro_func(lim_freq) + f_pro2_func(lim_freq))/20.0)
gain_correction = list(gain_correction) + list(gain_correction[::-1]) #反転して合わせる

fig, axes = plt.subplots(figsize=(8,4))
axes.plot(freq_array/1000.0, gain_correction)
axes.set_xlabel("Frequency kHz")
axes.set_ylabel("Relative response")
axes.set_ylim(-0.1, 1.1)

3. FIRフィルタの係数を計算する

波形データにフィルタ係数を畳み込み計算した結果は、それぞれのデータのフーリエ変換の積になります。フィルタ係数のフーリエ変換は、実現したい周波数特性そのものです。従って、先程の周波数特性を逆フーリエ変換すればフィルタの係数が求まります。また、FIRフィルタの係数は左右対称になります。

b = list(np.real(np.fft.ifft(gain_correction)))[0:N//2]
b = b[::-1] + b[1:]
b = b * np.hamming(N-1)
b /= np.sum(b) # 係数の和を1にする

フィルタ係数をフーリエ変換してみます。

4. FIRフィルタの適用

計算したFIRフィルタをホワイトノイズに適用して、適用された波形のスペクトルを確認してみます。

NN = 8192
ave_num = 1000
y = np.random.rand(NN*ave_num) # ホワイトノイズ
yy = signal.lfilter(b, 1, y) # FIRフィルタを適用した信号

fft_sum_w = np.zeros(NN)
fft_sum_f = np.zeros(NN)
for i in range(ave_num):
    fft_sum_w += np.abs(np.fft.fft(y[NN*i:NN*(i+1)]) / NN)
    fft_sum_f += np.abs(np.fft.fft(yy[NN*i:NN*(i+1)]) / NN)
fft_ave_w = 20.0 * np.log10(fft_sum_w / ave_num)
fft_ave_f = 20.0 * np.log10(fft_sum_f / ave_num)

fig, axes = plt.subplots(figsize=(8,4))
axes.plot(np.linspace(0.0, 500.0, num=NN), fft_ave_w)
axes.plot(np.linspace(0.0, 500.0, num=NN), fft_ave_f)
axes.set_xlim(1.0, 250.0)
axes.set_ylim(-50.0-15, -50.0)
axes.set_xscale("log")
axes.set_xlabel("Frequency kHz")
axes.set_ylabel("Amplitude dB")

ちゃんと実現したかった周波数特性になっています。

まとめ

実現したい周波数特性からFIRフィルタの係数を計算し、波形に適用しました。
今後はリアルタイム信号処理にも挑戦してみたいです。