n dB/octのスロープを持つフィルタ

n dB/oct のスロープを持つフィルタを設計します。

Image of the filter with n dB/oct slope.
A(ω)={(h/l)M,(πωh)(ω/l)M,(h<ω<l)1,(lωl)(ω/l)M,(l<ω<h)(h/l)M,(hωπ)M=m20log210 \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 &= {{\,m}\over{20}}\log_2 10 \end{aligned}

逆離散フーリエ変換を使う方法

逆離散フーリエ変換(IDFT)を使って楽に設計します。

コードを順番にPython3のインタープリタにコピペしていけば動作します。別ページでまとめて見ることもできます。

SciPyNumPyMatplotlibを 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の特性を作りました。

Image of filter specification.

フィルタ係数を出します。

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)

フィルタ係数の見た目です。

Image of the filter coefficients.

フィルタの周波数応答です。作った特性と一致します。

Image of frequency response of the filter.

不連続点を拡大するとギブス現象が観察できます。

showResponse(sample_rate, filt)
Image of gibbs phenomenon.

図の左は周波数軸が対数でないときのフィルタの周波数応答です。右上は低域側の不連続点、右下は高域側の不連続点を拡大したものです。

フィルタができたので、ノイズを作ってかけてみます。

# フィルタをノイズに適用。 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)

フィルタをかける前とかけた後のスペクトラムです。

Image of frequency response of the filter.

式を途中まで解く

解析解について調べます。

問題を再掲します。

A(ω)={(h/l)M,(πωh)(ω/l)M,(h<ω<l)1,(lωl)(ω/l)M,(l<ω<h)(h/l)M,(hωπ)M=m20log210 \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 &= {{\,m}\over{20}}\log_2 10 \end{aligned}

A(ω)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);
12πππA(ω)ejωndω=hMlMπheinωdω+1lMhlωMeinωdω+lleinωdω+1lMlhωMeinωdω+hMlMhπeinωdω=sin(ln)nπ+hMsin(nπ)lMnπhMsin(hn)lMnπ+i(1)M2nM12πlM(Γ(M+1,ihn)Γ(M+1,iln))+i(i)MnM12πlM(Γ(M+1,iln)Γ(M+1,ihn)) \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}

Γ(s,x)\Gamma(s, x)upper incomplete gamma function です。

Γ(s,x)=xts1etdt \Gamma(s, x) = \int_{x}^{\infty} t^{s - 1} e^{-t} dt

Γ(s,x)\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);
ππAP(ω)ejωndω=i(1)M2nM12πlM(Γ(M+1,ihn)Γ(M+1,iln))ππAR(ω)ejωndω=i(i)MnM12πlM(Γ(M+1,iln)Γ(M+1,ihn)) \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}

適当に Γ(M+1,ihn)\Gamma\left(M+1 , i\,h\,n\right) を展開します。

Γ(M+1,ihn)=ihntMetdt \Gamma\left(M+1 , i\,h\,n\right) = \int_{i\,h\,n}^{\infty} t^{M} e^{-t} dt

複素数の線積分に見えます。

参考サイト