one_pole_lowpass

Exponential moving average (EMA) フィルタではナイキスト周波数の成分が残るので、バイリニア変換した 1 次ローパスフィルタを作ることにしました。

以降では \(T=2\) とします。 \(T=2\) とするとカットオフ周波数を \(f_c / f_s\) と正規化できます。 \(f_c\) はカットオフ周波数、 \(f_s\) はサンプリング周波数です。

バイリニア変換の式を \(H(s)\) に代入します。また、簡略化のために \(k=1/\omega_c\) と置きます。

\[ \begin{aligned} H(z) &= \frac{1}{1 + k (1-z^{-1})/(1+z^{-1})} \\ &= \frac{1 + z^{-1}}{1+z^{-1} + k - k z^{-1}} \\ &= \frac{1 + z^{-1}}{(1 + k) + (1 - k) z^{-1}} \\ \end{aligned} \]

以下はバイリニア変換による周波数歪みを考慮した \(k\) の計算式です。

Python 3 で実装します。 sos は scipy.signal で使われる 2 次セクション (second order sections) の形式のフィルタ係数です。

import numpy as np

sampleRate = 48000
cutoffHz = 10000

k = 1 / np.tan(np.pi * cutoffHz / sampleRate)
a0 = 1 + k
a1 = (1 - k) / a0
b = 1 / a0
sos = [[b, b, 0, 1, a1, 0]]

以下は周波数特性のプロットです。ゲイン特性とカットオフ周波数が -3 dB で交差していることが確認できました。厳密な交差点のゲインは \(20 \log_{10} (1/\sqrt{2}) \approx -3.0103\ \text{dB}\) です。

ナイキスト周波数のゲインは \(z=-1\) を \(H(z)\) に代入すると得られます。 \(H(z)\) は分母が \(1 + z^{-1}\) なので理論上のゲインは 0 です。

実際に計算すると理論通りに動かないことがあるので以下のコードで確認しました。 x が入力、 y が出力です

x = np.ones(2**16)
x[1::2] = -1
y = signal.sosfilt(sos, x)

インパルス応答のような波形になっているのは、何もないところからサイン波がいきなり立ち上がったからです。 x が緩やかにフェードインするように窓関数をかけて再度試しました。

x = np.ones(2**16)
x[1::2] = -1
x *= signal.get_window("hann", len(x))  # 追加。
y = signal.sosfilt(sos, x)

以下は y のプロットです。図を拡大すると振幅は小さいものの完全には 0 にならないことが分かります。それでも EMA より良い性能を維持しています。

#include <cmath>
#include <numbers>

// BLT: Bilinear transform.
template<typename Sample> class BLTLP1 {
private:
  Sample bn = 1;
  Sample a1 = -1; // Negated.
  Sample x1 = 0;
  Sample y1 = 0;

public:
  void reset()
  {
    x1 = 0;
    y1 = 0;
  }

  void setCutoff(Sample sampleRate, Sample cutoffHz)
  {
    constexpr auto pi = std::numbers::pi_v<Sample>;
    auto k = Sample(1) / std::tan(pi * cutoffHz / sampleRate);
    auto a0 = Sample(1) + k;
    bn = Sample(1) / a0;
    a1 = (k - Sample(1)) / a0; // Negated.
  }

  Sample process(Sample x0)
  {
    auto y0 = bn * (x0 + x1) + a1 * y1;
    x1 = x0;
    return y1 = y0;
  }
};

\[ \begin{aligned} H(z) &= \frac{k \frac{1-z^{-1}}{1+z^{-1}}}{1 + k \frac{1-z^{-1}}{1+z^{-1}}} \\ &= \frac{k - k z^{-1}}{1 + z^{-1} + k - k z^{-1}} \\ &= \frac{k - k z^{-1}}{(1 + k) + (1 - k) z^{-1}} \\ \end{aligned} \]

k = 1 / np.tan(np.pi * cutoffHz / sampleRate)
a0 = 1 + k
b0 = k / a0
b1 = -b0
a1 = (1 - k) / a0
sos = [[b0, b1, 0, 1, a1, 0]]

以下はバイリニア変換による 1 次ハイパスの周波数応答です。

template<typename Sample> class BLTHP1 {
private:
  Sample b0 = 1;
  Sample a1 = 1;
  Sample x1 = 0;
  Sample y1 = 0;

public:
  void reset()
  {
    x1 = 0;
    y1 = 0;
  }

  void setCutoff(Sample sampleRate, Sample cutoffHz)
  {
    constexpr auto pi = std::numbers::pi_v<Sample>;
    auto k = Sample(1) / std::tan(pi * cutoffHz / sampleRate);
    auto a0 = Sample(1) + k;
    b0 = k / a0;
    a1 = (Sample(1) - k) / a0;
  }

  Sample process(Sample x0)
  {
    auto y0 = b0 * (x0 - x1) - a1 * y1;
    x1 = x0;
    return y1 = y0;
  }
};

\[ \begin{aligned} H(z) &= \frac{k_p}{1 + (k_p - 1) z^{-1}}, \\ k_p &= -y + \sqrt{y^2 + 2 y}, \\ y &= 1 - \cos(2 \pi f_c / f_s). \end{aligned} \]

以下は scipy.signal の sos 形式で EMA フィルタを設計するコードです。

def getEMASos(sampleRate, cutoffHz):
    y = 1 - np.cos(2 * np.pi * cutoffHz / sampleRate)
    kp = np.sqrt(y * y + 2 * y) - y
    return [[kp, 0, 0, 1, kp - 1, 0]]

\[ \begin{aligned} H_{\mathrm{HP}}(z) &= 1 - H(z) \\ &= 1 - \frac{k_p}{1 + (k_p - 1) z^{-1}} \\ &= \frac{(1 - k_p) + (k_p - 1) z^{-1}}{1 + (k_p - 1) z^{-1}} \\ \end{aligned} \]

def getEMASos(sampleRate, cutoffHz):
    y = 1 - np.cos(2 * np.pi * cutoffHz / sampleRate)
    kp = np.sqrt(y * y + 2 * y) - y
    kq = 1 - kp
    return [[kq, -kq, 0, 1, -kq, 0]]