お手軽なFIRフィルタのレシピ

FIRのローパス、ハイパス、バンドパス、バンドリジェクトフィルタを作ります。

記号

通過域の大きさ(Amplitude)は1とします。

ローパスフィルタ

Image of FIR lowpass filter.

A(ω)={1,(ωlωωl)0,otherwise A(\omega) = \begin{cases} 1, & (-\omega_l \leq \omega \leq \omega_l) \\ 0, & \text{otherwise} \end{cases}

A(ω)A(\omega)逆離散時間フーリエ変換の式に入れて解きます。

12πππA(ω)ejωndω=12πωlωlejωndω=sin(ωln)nπ \begin{aligned} \frac{1}{2\pi}\int^{\pi}_{-\pi} A(\omega) e^{j\omega n} d\omega &= \frac{1}{2\pi}\int^{\omega_l}_{-\omega_l} e^{j\omega n} d\omega \\ &= \frac{\sin (\omega_l n)}{n \pi} \end{aligned}

コードに変えます。

function makeLowpassCoefficient(length, cutoff) { // cutoff の範囲は [0, 1] var coefficient = new Array(length).fill(0) var half = (length % 2 === 0 ? length - 1 : length) / 2 var omegaL = 2 * Math.PI * cutoff for (var i = 0; i < length; ++i) { var n = i - half coefficient[i] = (n === 0) ? 1 : Math.sin(omegaL * n) / (Math.PI * n) } return coefficient }

ハイパスフィルタ

Image of FIR highpass filter.

A(ω)={1,(ωsωωh)1,(ωhωωs)0,otherwise A(\omega) = \begin{cases} 1, & (-\omega_s \leq \omega \leq -\omega_h) \\ 1, & (\omega_h \leq \omega \leq \omega_s) \\ 0, & \text{otherwise} \end{cases}

A(ω)A(\omega) を逆離散時間フーリエ変換の式に入れて解きます。

12πππA(ω)ejωndω=12π(πωhejωndω+ωhπejωndω)=sin(nπ)nπsin(ωhn)nπ \begin{aligned} \frac{1}{2\pi}\int^{\pi}_{-\pi} A(\omega) e^{j\omega n} d\omega &= \frac{1}{2\pi} \biggl ( \int^{-\omega_h}_{-\pi} e^{j\omega n} d\omega + \int^{\pi}_{\omega_h} e^{j\omega n} d\omega \biggr ) \\ &= {{\sin \left(n\,\pi\right)}\over{n\,\pi}}-{{\sin \left(\omega_h\,n\right)}\over{n\,\pi}} \end{aligned}

コードに変えます。

function makeHighpassCoefficient(length, cutoff) { // cutoff の範囲は [0, 1] var coefficient = new Array(length).fill(0) var half = (length % 2 === 0 ? length - 1 : length) / 2 var omegaH = 2 * Math.PI * cutoff for (var i = 0; i < length; ++i) { var n = i - half coefficient[i] = (n === 0) ? 1 : (Math.sin(Math.PI * n) - Math.sin(omegaH * n)) / (Math.PI * n) return coefficient }

バンドパスフィルタ

Image of FIR bandpass filter.

A(ω)={1,(ωhωωl)1,(ωlωωh)0,otherwise A(\omega) = \begin{cases} 1, & (-\omega_h \leq \omega \leq -\omega_l) \\ 1, & (\omega_l \leq \omega \leq \omega_h) \\ 0, & \text{otherwise} \end{cases}

A(ω)A(\omega) を逆離散時間フーリエ変換の式に入れて解きます。

12πππA(ω)ejωndω=12π(ωhωlejωndω+ωlωhejωndω)=sin(ωhn)nπsin(ωln)nπ \begin{aligned} \frac{1}{2\pi}\int^{\pi}_{-\pi} A(\omega) e^{j\omega n} d\omega &= \frac{1}{2\pi} \biggl ( \int^{-\omega_l}_{-\omega_h} e^{j\omega n} d\omega + \int^{\omega_h}_{\omega_l} e^{j\omega n} d\omega \biggr ) \\ &= {{\sin \left(\omega_h\,n\right)}\over{n\,\pi}}-{{\sin \left(\omega_l\,n\right)}\over{n\,\pi}} \end{aligned}

コードに変えます。

function makeBandpassCoefficient(length, low, high) { // low, high の範囲は [0, 1] var coefficient = new Array(length).fill(0) var half = (length % 2 === 0 ? length - 1 : length) / 2 var twoPi = 2 * Math.PI var omegaL = twoPi * low var omegaH = twoPi * high for (var i = 0; i < length; ++i) { var n = i - half coefficient[i] = (n === 0) ? 1 : (Math.sin(omegaH * n) - Math.sin(omegaL * n)) / (Math.PI * n) return coefficient }

バンドリジェクトフィルタ

バンドストップフィルタとも呼ばれるようです。

Image of FIR bandreject filter.

A(ω)={1,(ωsωωh)1,(ωlωωl)1,(ωhωωs)0,otherwise A(\omega) = \begin{cases} 1, & (-\omega_s \leq \omega \leq -\omega_h) \\ 1, & (-\omega_l \leq \omega \leq \omega_l) \\ 1, & (\omega_h \leq \omega \leq \omega_s) \\ 0, & \text{otherwise} \end{cases}

A(ω)A(\omega) を逆離散時間フーリエ変換の式に入れて解きます。

12πππA(ω)ejωndω=12π(πωhejωndω+ωlωlejωndω+ωhπejωndω)=sin(nπ)nπ+sin(ωln)nπsin(ωhn)nπ \begin{aligned} \frac{1}{2\pi}\int^{\pi}_{-\pi} A(\omega) e^{j\omega n} d\omega &= \frac{1}{2\pi} \biggl ( \int^{-\omega_h}_{-\pi} e^{j\omega n} d\omega + \int^{\omega_l}_{-\omega_l} e^{j\omega n} d\omega + \int^{\pi}_{\omega_h} e^{j\omega n} d\omega \biggr ) \\ &= {{\sin \left(n\,\pi\right)}\over{n\,\pi}}+{{\sin \left(\omega_l\,n\right)}\over{n\,\pi}}-{{\sin \left(\omega_h\,n\right)}\over{n\,\pi}} \end{aligned}

コードに変えます。

function makeBandrejectCoefficient(length, low, high) { // low, high の範囲は [0, 1] var coefficient = new Array(length).fill(0) var half = (length % 2 === 0 ? length - 1 : length) / 2 var twoPi = 2 * Math.PI var omegaL = twoPi * low var omegaH = twoPi * high for (var i = 0; i < length; ++i) { var n = i - half var piN = Math.PI * n coefficient[i] = (n === 0) ? 1 : (Math.sin(piN) + Math.sin(omegaL * n) - Math.sin(omegaH * n)) / piN return coefficient }

FIRフィルタのかけ方

まずフィルタをかけるソースを用意します。

var audioContext = new AudioContext() var source = new Array(audioContext.sampleRate) for (var i = 0; i < source.length; ++i) { source[i] = 2.0 * Math.random() - 1.0 // 適当なノイズ。 }

FIRフィルタは適切な窓関数をかけることで特性が改善します。ここではお手軽でそれなりに特性がいいBlackman–Harris窓を使います。

function blackmanHarrisWindow(length) { var window = new Array(length) var a0 = 0.35875 var a1 = 0.48829 var a2 = 0.14128 var a3 = 0.01168 var pi_N1 = Math.PI / (window.length - 1) var twopi_N1 = 2 * pi_N1 var fourpi_N1 = 4 * pi_N1 var sixpi_N1 = 6 * pi_N1 for (var n = 0; n < window.length; ++n) { window[n] = a0 - a1 * Math.cos(n * twopi_N1) + a2 * Math.cos(n * fourpi_N1) - a3 * Math.cos(n * sixpi_N1) } return window }

フィルタを用意します。

var filterLength = 1025 var cutoff = 1000 // Hz var coefficient = makeLowpassCoefficient(filterLength, cutoff / audioContext.sampleRate) var filter = blackmanHarrisWindow(filterLength).map((v, i) => v * coefficient[i])

畳み込み(Convolution)を行ってソースにFIRフィルタをかけます。

var destination = new Array(source.length).fill(0) var buffer = new Array(filter.length).fill(0) for (var i = 0; i < source.length; ++i) { buffer.push(source[i]) buffer.shift() for (var j = 0; j < filter.length; ++j) { destination[i] += buffer[j] * filter[j] } }

Computer Algebra System の利用

手で式を解くと間違えることがあるので Computer Algebra System (CAS) を利用します。

今回のような簡単な式であればWolfram Alphaが便利です。Wolfram Alphaでは数字でない下付き文字が使えないようなので ωl\omega_lll に置き換えています。以降のCASのコードも同じ置き換えを使います。

(integral e^(i*omega*n) for omega from -l to l) / (2pi)

Maximaは式の整理について指定する必要があります。コードの demoivreexp(%i*n)%i * sin(n) + cos(n) に置き換えています。

expand(demoivre(integrate(exp(%i * omega * n) / (2 * pi), omega, -l, l)));

SymPyも使えますが少し長めです。 rewrite(sin) でオイラーの公式を適用しています。 n = Symbol('n', positive=True) が無いとコードの integrate(...) を解いてくれません。

# SymPy 1.1.1 from sympy import * n = Symbol('n', positive=True) l = Symbol('l') omega = Symbol('omega') answer = simplify(integrate(exp(I * omega * n), (omega, -l, l)).rewrite(sin)) pprint(answer)

フィルタ係数の計算について

フィルタ係数は周波数特性を逆離散時間フーリエ変換することで得られます。

12πππA(ω)ejωndω \frac{1}{2\pi}\int^{\pi}_{-\pi} A(\omega) e^{j\omega n} d\omega

この逆離散時間フーリエ変換の式に、周波数ffを対応させて使うときはπf/fs\pi{f}/{f_s}と変換します。

逆離散時間フーリエ変換をするときは[fs, fs][-f_s,\ f_s]の範囲を考慮する必要がありますが、ここまでに出てきた周波数特性の図では[fs, 0][-f_s,\ 0]の範囲を省略していました。[fs, 0][-f_s,\ 0]の範囲での周波数特性は[0, fs][0,\ f_s]の鏡像になっています。

ローパスフィルタを例に見ていきます。

Image of range [-f_s, f_s] of lowpass frequency responce.

偶関数の積分の性質を利用して式を変形できそうなので試してみます。

12πππA(ω)ejωndω=12πωlωlejωndω=22π0ωlejωndω=sin(ωln)nπicos(ωln)nπ+inπ \begin{aligned} \frac{1}{2\pi}\int^{\pi}_{-\pi} A(\omega) e^{j\omega n} d\omega &= \frac{1}{2\pi}\int^{\omega_l}_{-\omega_l} e^{j\omega n} d\omega \\ &= \frac{2}{2\pi}\int^{\omega_l}_{0} e^{j\omega n} d\omega \\ &= {{\sin \left(\omega_l\,n\right)}\over{n\,\pi}}-{{i\,\cos \left(\omega_l\,n\right)}\over{n\,\pi}}+{{i}\over{n\,\pi}} \end{aligned}

虚部が出てきました。実部は式を変形をしないときと同じになっています。

念のために定義どおり計算した方がよさそうです。

その他

式を解かなくても周波数特性を逆離散フーリエ変換すればフィルタは作れます。

フィルタ係数が固定のときはSciPyOctaveなどを使って設計するほうが楽で確実です。

ここで作ったバンドパスフィルタを Banded Waveguides に使おうとしたのですが、切れ味が良すぎてWaveguide間の干渉がほとんど起こらず、面白い音になりませんでした。

参考サイト