この記事では、
畳み込み演算の実装
畳み込みの定義
入力信号と、システムの因果性を考慮して以下の式で実装します。
$y(n)= \sum^{n-1}_{k=0}{x(k)h(n-k)}$
これ元に参考になるべく繰り返しの回数が減るようにすると
def conv1d(x, h):
"""
convolve x and h
Parameters
----------
x : array_like
First one-dimensional input array.
h : array_like
Second one-dimensional input array.
Returns
-------
out : ndarray
Discrete, linear convolution of a and v.
"""
N = len(x)
out = np.zeros(N + len(h))
for n in range(len(h)):
out[n: n + N] += x * h[n]
return out[:N]
こんな感じになります。
ディジタルフィルタの設計
バンドストップフィルタの設計
基本的にどんなフィルタでも考え方は同じです。
今回はバンドストップフィルタで行います。
Band Stop Filter を設計するためには周波数領域で、以下のようになるインパルス応答を実装すればよいです。(ただし$ω$は正規化された角周波数)



先の式を離散逆フーリエ変換した結果は以下のようになります。(ただし、 nは離散信号の1点を表す)



これがバンドストップフィルタのインパルス応答です。
なので、必要な分の窓幅$N$を定義して、そのサイズのバンドストップフィルタを返す関数を作成していきます。
バンドストップフィルタのインパルス応答
自作sinc関数を定義
numpyのsinc関数では定義が
$\frac{1}{πx}sin(πx)$となっており、式変形がわずらわしいので、自分で定義します。
ただ、自分で定義したことにより、入力が要素しか受け付けなくなっている点に注意が必要です。
def sinc(x):
"""
sinc function
Parameters
----------
x : float
.input
Returns
-------
y : float
apply sinc function .
"""
if x == 0.0:
return 1.0
else:
return np.sin(x) / x
入力のインデックスに対して、返答するインパルス応答を定義
from math import pi
def BSF_impulse(index, f_low, f_high, sr):
"""
calculate BSF impulse at each point
Parameters
----------
idx : int
index number
f_low : float
low edge frequency
f_high : float
high edge frequency
sr : int
sampling rate
Returns
-------
y : float
BSF impulse
"""
# normalize sampling rate
f_low = 2 * pi * f_low / sr
f_high = 2 * pi * f_high / sr
return sinc(pi * index) + (f_low * sinc(f_low * index) - f_high * sinc(f_high * index)) / pi
バンドストップフィルタを作成
今回は、窓関数法で適応をします。
フィルタは、0を中心に奇関数として定義をされているので、そうなるようにインデックスを入力します。
しかし、実際に使う際は、$(0, N+1)$の範囲で定義されていた方が使いやすいので、
b.append()のように逐次追加をしていき、最後にnumpy配列に変換をしています。
from scipy import hamming
def BandStopFilter(f_low, f_high, sr, N=30):
"""
make band stop filter
Parameters
----------
f_low : float
low edge frequency
f_high : float
high edge frequency
N : int, default 30
the number of filter coefficients
sr : float
sampling rate
Returns
-------
b : ndarray
Filter coefficient
"""
if N % 2 != 0:
N += 1
N = int(N)
b = []
# calculate filter coefficients
for idx in range(-N // 2, N // 2 + 1):
b.append(BSF_impulse(idx, f_low, f_high, sr))
# Apply window function method
b = np.array(b)
window_func = hamming(N + 1)
return b * window_func
上で定義したBFS_impulse()関数を変更すれば任意のフィルタを設計することができます。
今回作成したフィルタは、
窓幅:$N = 100$
阻止域周波数:$f_{low} = 8000, f_{high} = 12000$
として定義をしました。
インパルス応答を確認



こんなフィルタが作成できます。
位相特性を確認



振幅特性を確認



どちらも定義した阻止域で周波数が遮断されていることが分かります。
実際に関数を適応
下の関数では、バンドストップフィルタ適応前後の周波数特性と位相特性と、バンドストップフィルタ適応前後のスペクトログラムを表示しています。
import matplotlib.pyplot as plt
import librosa
def main():
N = args.N # [step]
f_low = args.f_low # [Hz]
f_high = args.f_high # [Hz]
wav, sr = librosa.load(fname, sr=48000) # [sec], [Hz]
x = wav
F_x = np.fft.fft(x * hamming(len(x)))
amp = 20 * np.log10(np.abs(F_x))
arg = np.unwrap(np.angle(F_x)) * 180 / pi
frequency_range = np.arange(0, sr // 2, (sr // 2) / (wav.size // 2 + 1))
plt.figure(figsize=(15, 5))
plt.subplot(1, 2, 1)
plt.title('Amplitude characteristic', size=16)
plt.plot(frequency_range, amp[: len(amp) // 2 + 1])
plt.xlabel("Frequency [Hz]", size=14)
plt.ylabel("Amplitude [dB]", size=14)
plt.subplot(1, 2, 2)
plt.title('Phase characteristic', size=16)
plt.plot(frequency_range, arg[: len(arg) // 2 + 1])
plt.xlabel("Frequency [Hz]", size=14)
plt.ylabel("Phase [rad]", size=14)
plt.savefig("Original_Amp_Phase_char.png")
plt.close()
b = BandStopFilter(f_low, f_high, sr, N=N)
x_hat = conv1d(x, b)
F_x = np.fft.fft(x_hat * hamming(len(x_hat)))
amp = 20 * np.log10(np.abs(F_x))
arg = np.unwrap(np.angle(F_x)) * 180 / pi
frequency_range = np.arange(0, sr // 2, (sr // 2) / (wav.size // 2 + 1))
plt.figure(figsize=(15, 5))
plt.subplot(1, 2, 1)
plt.title('Amplitude characteristic', size=16)
plt.plot(frequency_range, amp[: len(amp) // 2 + 1])
plt.xlabel("Frequency [Hz]", size=14)
plt.ylabel("Amplitude [dB]", size=14)
plt.subplot(1, 2, 2)
plt.title('Phase characteristic', size=16)
plt.plot(frequency_range, arg[: len(arg) // 2 + 1])
plt.xlabel("Frequency [Hz]", size=14)
plt.ylabel("Phase [rad]", size=14)
plt.savefig("BST_Amp_Phase_char.png")
plt.close()
plt.figure(figsize=(12, 5))
plt.title("Original signal spectrogram", size=16)
plt.xlabel("Time [sec]", size=14)
plt.ylabel("Frequency [Hz]", size=14)
spectrum, freq, _, _ = plt.specgram(x, Fs=sr)
cbar = plt.colorbar()
cbar.set_label("magnitude [dB]", size=14)
plt.savefig("Original_signal_spectrogram.png")
plt.close()
plt.figure(figsize=(12, 5))
plt.title("Band Stop Filter signal spectrogram", size=16)
plt.xlabel("Time [sec]", size=14)
plt.ylabel("Frequency [Hz]", size=14)
_=plt.specgram(x_hat, Fs=sr)
cbar = plt.colorbar()
cbar.set_label("magnitude [dB]", size=14)
plt.savefig("BSF_signal_spectrogram.png")
plt.close()
librosa.output.write_wav('BSF_signal.wav', x_hat, sr)
print("over")
if __name__ == "__main__":
main()
元の波の振幅特性と位相特性



BSF適応後の振幅特性と位相特性



位相特性はわかりにくいですが、振幅特性は削られていることが分かります。
元の波のスペクトログラム



BSF後のスペクトログラム



阻止域周波数:$f_{low} = 8000, f_{high} = 12000$の部分がごっそり取り除かれていることが分かります。
目的のフィルタが完成したと考えて良いのではないでしょうか!
最後に
今回はバンドストップフィルタでしたが、同じ方法でどんなフィルタでも作成できることが分かりました。
手元でフィルタを作成しているとCNNがどういった働きをしているのかがよく分かります。
すごく荒く説明をすると、目的の周波数を強調するフィルタを作成して、その周波数成分をもとに回帰問題を解いているのです。
何事も基礎から学んでいくことは大切だとしみじみ感じました!






タイピングもままならない完全にプログラミング初心者から
たった二ヶ月で
応用も簡単にできる…!!
という状態になるまで、一気に成長させてくれたオススメのプログラミングスクールをご紹介します!


