何かあれば GitHub のリポジトリに issue を作るか ryukau@gmail.com までお気軽にどうぞ。
Update: 2023-06-22
n dB/oct のスロープを持つフィルタを設計します。
\[ \begin{aligned} A(\omega) &= \begin{cases} (h/l)^{M}, & (-\pi \leq \omega \leq -h) \\ (\omega/l)^{M}, & (-h \lt \omega \lt -l) \\ 1, & (-l \leq \omega \leq l) \\ (\omega/l)^{M}, & (l \lt \omega \lt h) \\ (h/l)^{M}, & (h \leq \omega \leq \pi) \end{cases} \\\\ M &= {{\,n}\over{20}}\log_2 10 \end{aligned} \]
逆離散フーリエ変換（IDFT）を使って楽に設計します。
コードを順番にPython3のインタープリタにコピペしていけば動作します。別ページでまとめて見ることもできます。
SciPy、NumPy、Matplotlibを import します。
import matplotlib.pyplot as p
import numpy as np
from scipy import signal
from scipy.fftpack import *プロットに使う関数を用意します。
def showResponse(sample_rate, sig):
# sig の周波数成分を dB に変換してプロット。
res = fftshift(fft(sig, 2**16) / (len(sig) / 2.0))
gain = np.array_split(
20 * np.log10(abs(res / np.max(abs(res)))),
2,
)[0][::-1]
freq = getFreq(sample_rate, gain)
p.plot(freq, gain)
p.grid(which='both')
p.show()
def getFreq(sample_rate, spec):
# 出力はプロットのx軸に使う。
freq_step = sample_rate / 2 / (len(spec) - 1)
return [i * freq_step for i in range(0, len(spec))]
def toDecibel(spec):
return [20 * np.log10(v) for v in spec]欲しい周波数特性を用意します。
def makeFrequencySpecification(sample_rate, length, low, high, slope):
# slope は負の値。単位は[dB/oct]。
# low, high は周波数[Hz]。
#
# dB/oct を算出する式の解説。
# https://femci.gsfc.nasa.gov/random/randomequations.html
#
# 6.020599913279624 = 20 * np.log10(2)
#
oct = (high / low)**(slope / 6.020599913279624)
spec = [0] * length
freq_step = sample_rate / 2 / (length - 1)
for i in range(0, len(spec)):
freq = i * freq_step
if freq <= low:
spec[i] = 1
elif freq >= high:
spec[i] = oct
else:
spec[i] = (freq / low)**(slope / 6.020599913279624)
return spec
sample_rate = 44100
spec = makeFrequencySpecification(sample_rate, 1024, 100, 10000, -10)
# 設計した特性をプロット。
spec_db = toDecibel(spec)
freq = getFreq(sample_rate, spec)
p.title('Specification')
p.plot(freq, spec_db)
p.xscale('log')
p.grid(which='both')
p.show()ここでは100Hzから10000Hz間において、-10dB/octの傾きで減衰する、窓長1024の特性を作りました。
フィルタ係数を出します。
filt = fftshift(ifft(spec + spec[::-1]))
# フィルタ係数の見た目。
p.title('Filter Coefficients')
p.plot(filt)
p.show()
# フィルタの周波数応答。
p.title('Frequency Response of the filter')
p.xscale('log')
showResponse(sample_rate, filt)フィルタ係数の見た目です。
フィルタの周波数応答です。作った特性と一致します。
不連続点を拡大するとギブス現象が観察できます。
showResponse(sample_rate, filt)
図の左は周波数軸が対数でないときのフィルタの周波数応答です。右上は低域側の不連続点、右下は高域側の不連続点を拡大したものです。
フィルタができたので、ノイズを作ってかけてみます。
# フィルタをノイズに適用。
noise = np.random.uniform(-1.0, 1.0, sample_rate)
p.title('Before')
p.xscale('log')
showResponse(sample_rate, noise)
filtered = signal.convolve(noise, filt, mode='same')
p.title('After')
p.xscale('log')
showResponse(sample_rate, filtered)フィルタをかける前とかけた後のスペクトラムです。
解析解について調べます。
問題を再掲します。
\[ \begin{aligned} A(\omega) &= \begin{cases} (h/l)^{M}, & (-\pi \leq \omega \leq -h) \\ (\omega/l)^{M}, & (-h \lt \omega \lt -l) \\ 1, & (-l \leq \omega \leq l) \\ (\omega/l)^{M}, & (l \lt \omega \lt h) \\ (h/l)^{M}, & (h \leq \omega \leq \pi) \end{cases} \\\\ M &= {{\,n}\over{20}}\log_2 10 \end{aligned} \]
\(A(\omega)\) を逆DTFTの式に入れてMaximaで解きます。
/* Maxima */
expintrep: true;
expintexpand: true;
declare(n, integer, M, real, l, real, h, real);
assume(n > 0, l >= 0, l <= pi, h >= 0, h <= pi, h > l);
oct(m, l, h) := (h / l)**(M);
O: integrate(oct(m, l, h) * exp(%i*omega*n), omega, -pi, -h);
P: integrate(oct(m, l, omega) * exp(%i*omega*n), omega, -h, -l);
Q: integrate(exp(%i*omega*n), omega, -l, l);
R: integrate(oct(m, l, omega) * exp(%i*omega*n), omega, l, h);
S: integrate(oct(m, l, h) * exp(%i*omega*n), omega, h, pi);
answer: expand(demoivre((O + P + Q + R + S) / (2*pi)));
/* where */
M: log(10)*m/20/log(2);\[ \begin{aligned} \frac{1}{2\pi}\int^{\pi}_{-\pi} A(\omega) e^{j\omega n} d\omega = &{{h^{M}}\over{l^{M}}}\int_{-\pi}^{-h}{e^{i\,n\,\omega}\;d\omega}\\ &+{{1}\over{l^{M}}}\int_{-h}^{-l}{\omega^{M}\,e^{i\,n\,\omega}\;d\omega}\\ &+\int_{-l}^{l}{e^{i\,n\,\omega}\;d\omega}\\ &+{{1}\over{l^{M}}}\int_{l}^{h}{\omega^{M}\,e^{i\,n\,\omega}\;d\omega}\\ &+{{h^{M}}\over{l^{M}}}\int_{h}^{\pi}{e^{i\,n\,\omega}\;d\omega}\\ = &{{\sin \left(l\,n\right)}\over{n\,\pi}} +{{h^{M}\,\sin \left(n\,\pi\right)}\over{l^{M}\,n\,\pi}} -{{h^{M}\,\sin \left(h\,n\right)}\over{l^{M}\,n\,\pi}} \\ &+{{i\,\left(-1\right)^{{{M}\over{2}}} n^{-M-1}}\over{2\,\pi\,l^{M}}} \Bigl(\Gamma\left(M+1 , i\,h\,n\right) - \Gamma\left(M+1 , i\,l\,n\right)\Bigr) \\ &+{{i\,\left(-i\right)^{-M}\,n^{-M-1}}\over{2\,\pi\,l^{M}}} \Bigl(\Gamma\left(M+1 , -i\,l\,n\right) - \Gamma\left(M+1 , -i\,h\,n\right)\Bigr) \end{aligned} \]
\(\Gamma(s, x)\) は upper incomplete gamma function です。
\[ \Gamma(s, x) = \int_{x}^{\infty} t^{s - 1} e^{-t} dt \]
\(\Gamma(s, x)\) がどこから出てきたのかを確認するために、Maximaのコードで定義したP式とR式に注目します。
oct(m, l, h) := (h / l)**(M);
assume(n > 0, l >= 0, l <= pi, h >= 0, h <= pi, h > l);
P: integrate(oct(m, l, omega) * exp(%i*omega*n), omega, -h, -l);
R: integrate(oct(m, l, omega) * exp(%i*omega*n), omega, l, h);\[ \begin{aligned} \int^{\pi}_{-\pi} A_P(\omega) e^{j\omega n} d\omega = &{{i\,\left(-1\right)^{{{M}\over{2}}} n^{-M-1}}\over{2\,\pi\,l^{M}}} \Bigl(\Gamma\left(M+1 , i\,h\,n\right) - \Gamma\left(M+1 , i\,l\,n\right)\Bigr) \\ \int^{\pi}_{-\pi} A_R(\omega) e^{j\omega n} d\omega = &{{i\,\left(-i\right)^{-M}\,n^{-M-1}}\over{2\,\pi\,l^{M}}} \Bigl(\Gamma\left(M+1 , -i\,l\,n\right) - \Gamma\left(M+1 , -i\,h\,n\right)\Bigr) \end{aligned} \]
適当に \(\Gamma\left(M+1 , i\,h\,n\right)\) を展開します。
\[ \Gamma\left(M+1 , i\,h\,n\right) = \int_{i\,h\,n}^{\infty} t^{M} e^{-t} dt \]
複素数の線積分に見えます。