ダウンサンプリング

import scipy.signal as signal
import numpy as np

def decimationFir(length, oversample, cutoff, nyquist=0.5):
    """オーバーラップありのときは nyquist を 0.6 などに設定する。"""
    bands = np.hstack((
        [0, cutoff],
        [nyquist, oversample / 2],
    ))
    desired = (1, 0)
    return signal.remez(length, bands, desired, fs=oversample)

fir = decimationFir(68, 4, 0.4, 0.6)

fir = decimationFir(108, 4, 0.375)

#include <array>
#include <iostream>
#include <numeric>
#include <vector>

/*
2 倍のオーバーサンプリングに特化した FIR フィルタ係数。

パスバンドに ±0.003 dB のリップル。
ストップバンドで -69 dB 低減。
オーバーラップ無し。

---python
import scipy.signal as signal
oversample = 2
signal.remez(64, [0, 0.375, 0.5, oversample], [1, 0], fs=oversample)
---
*/
template<typename Sample> struct FirRemez64Decimation2 {
  constexpr static size_t bufferSize = 32;
  constexpr static size_t nPhase = 2;

  constexpr static std::array<std::array<Sample, bufferSize>, nPhase> coefficient{{
    {Sample(-0.00012586358900962356), Sample(0.0001828145651312807),
     Sample(8.205416319885108e-05),   Sample(-0.0007239553549487794),
     Sample(0.0018923639640423555),   Sample(-0.00353192867622875),
     Sample(0.005331459582347676),    Sample(-0.006667992332855854),
     Sample(0.006636334343333397),    Sample(-0.004162169091786725),
     Sample(-0.0018313041529962014),  Sample(0.012291623826628002),
     Sample(-0.02812500771752753),    Sample(0.05102267479538475),
     Sample(-0.08758235952225661),    Sample(0.18605559261058696),
     Sample(0.4035431721330982),      Sample(-0.03627960335535102),
     Sample(-0.006631516329342658),   Sample(0.02064755543805059),
     Sample(-0.023817278019366648),   Sample(0.021451149873294335),
     Sample(-0.0164261830615856),     Sample(0.010706859341764241),
     Sample(-0.005602179907252967),   Sample(0.001817817136086864),
     Sample(0.00046580645271121403),  Sample(-0.0014510531292568874),
     Sample(0.0015411334474003102),   Sample(-0.0011736418939169283),
     Sample(0.0006849141024872853),   Sample(-0.00039292037784520776)},
    {Sample(-0.00039292037784520776), Sample(0.0006849141024872853),
     Sample(-0.0011736418939169283),  Sample(0.0015411334474003102),
     Sample(-0.0014510531292568874),  Sample(0.00046580645271121403),
     Sample(0.001817817136086864),    Sample(-0.005602179907252967),
     Sample(0.010706859341764241),    Sample(-0.0164261830615856),
     Sample(0.021451149873294335),    Sample(-0.023817278019366648),
     Sample(0.02064755543805059),     Sample(-0.006631516329342658),
     Sample(-0.03627960335535102),    Sample(0.4035431721330982),
     Sample(0.18605559261058696),     Sample(-0.08758235952225661),
     Sample(0.05102267479538475),     Sample(-0.02812500771752753),
     Sample(0.012291623826628002),    Sample(-0.0018313041529962014),
     Sample(-0.004162169091786725),   Sample(0.006636334343333397),
     Sample(-0.006667992332855854),   Sample(0.005331459582347676),
     Sample(-0.00353192867622875),    Sample(0.0018923639640423555),
     Sample(-0.0007239553549487794),  Sample(8.205416319885108e-05),
     Sample(0.0001828145651312807),   Sample(-0.00012586358900962356)},
  }};
};

template<typename Sample, typename FIR> class FirFilter {
  std::array<std::array<Sample, FIR::bufferSize>, FIR::nPhase> buf;

public:
  FirFilter() { reset(); }

  void reset()
  {
    for (auto &bf : buf) bf.fill(Sample(0));
  }

  Sample process(std::array<Sample, FIR::nPhase> &input)
  {
    Sample sum = 0;

    for (size_t ph = 0; ph < FIR::nPhase; ++ph) {
      auto &bf = buf[ph];
      for (size_t i = 0; i < bf.size() - 1; ++i) bf[i] = bf[i + 1];
      bf.back() = input[ph];

      const auto &fir = FIR::coefficient[ph];
      for (size_t i = 0; i < fir.size(); ++i) sum += bf[i] * fir[i];
    }

    return sum;
  }
};

int main()
{
  FirFilter<float, FirRemez64Decimation2<float>> filter;

  std::vector<float> data(64);
  std::iota(data.begin(), data.end(), 0.0f);

  std::array<float, 2> input;
  for (size_t i = 0; i < data.size(); i += 2) {
    input[0] = data[i];
    input[1] = data[i + 1];
    std::cout << filter.process(input) << "\n";
  }

  return 0;
}

auto &bf = buf[ph];
for (size_t i = 0; i < bf.size() - 1; ++i) bf[i] = bf[i + 1];
bf.back() = input[ph];

auto &bf = buf[ph];
std::rotate(buf.begin(), buf.begin() + 1, buf.end());
bf.back() = input[ph];

import scipy.signal as signal
import numpy as np

order = 12        # フィルタの次数。
ripple = 0.01     # パスバンドのリップルの大きさ。単位は dB 。
attenuation = 100 # ストップバンドでの低減。単位は dB 。
cutoff = 0.4      # カットオフ周波数。 fs = 1 に正規化された周波数。
oversample = 2    # オーバーサンプリングの倍率。

sos = signal.ellip(
    order, ripple, attenuation, cutoff, "low", output="sos", fs=oversample)

import scipy.signal as signal
import numpy as np

order = 8        # フィルタの次数。
cutoff = 0.4      # カットオフ周波数。 fs = 1 に正規化された周波数。
oversample = 2    # オーバーサンプリングの倍率。

sos = signal.butter(order, cutoff, "low", output="sos", fs=oversample)

#include <array>
#include <iostream>
#include <numeric>
#include <vector>

/*
---python
import numpy
from scipy import signal
sos = signal.ellip(12, 0.01, 100, 0.4, "low", output="sos", fs=8)
---
*/
template<typename Sample> struct SosEllipticDecimation8 {
  constexpr static size_t nSection = 6;

  constexpr static std::array<std::array<Sample, 5>, nSection> co{{
    {Sample(2.0478175697515062e-05), Sample(5.313899743869536e-06),
     Sample(2.047817569751506e-05), Sample(-1.7517115385026274),
     Sample(0.7700392502488301)},
    {Sample(1.0), Sample(-1.438980096034124), Sample(1.0), Sample(-1.774966834164871),
     Sample(0.8105614413580351)},
    {Sample(1.0), Sample(-1.7274028784750513), Sample(1.0), Sample(-1.808053369266683),
     Sample(0.867640704095363)},
    {Sample(1.0), Sample(-1.8120823742998813), Sample(1.0000000000000004),
     Sample(-1.8388054096377313), Sample(0.9192163691033324)},
    {Sample(1.0), Sample(-1.8441397779182482), Sample(1.0000000000000002),
     Sample(-1.8637392843186766), Sample(0.9580528037069385)},
    {Sample(1.0), Sample(-1.855844538281174), Sample(1.0), Sample(-1.8856454058448944),
     Sample(0.9871065355804314)},
  }};
};

// SOS: Second order sections.
template<typename Sample, typename IIR, size_t oversample = 8> class SosFilter {
public:
  void reset()
  {
    x1.fill(0);
    x2.fill(0);
    y1.fill(0);
    y2.fill(0);
  }

  void push(Sample input)
  {
    for (size_t i = 0; i < IIR::nSection; ++i) {
      // clang-format off
      auto y0 =
        + IIR::co[i][0] * input
        + IIR::co[i][1] * x1[i]
        + IIR::co[i][2] * x2[i]
        - IIR::co[i][3] * y1[i]
        - IIR::co[i][4] * y2[i];
      // clang-format on

      x2[i] = x1[i];
      x1[i] = input;
      y2[i] = y1[i];
      y1[i] = y0;
      input = y0;
    }
  }

  inline Sample output() { return y1[IIR::nSection - 1]; }

  Sample process(const std::array<Sample, oversample> &input)
  {
    for (const auto &value : input) push(value);
    return output();
  }

  std::array<Sample, IIR::nSection> x1{};
  std::array<Sample, IIR::nSection> x2{};
  std::array<Sample, IIR::nSection> y1{};
  std::array<Sample, IIR::nSection> y2{};
};

int main()
{
  SosFilter<float, SosEllipticDecimation8<float>> filter;

  std::vector<float> data(64);
  std::iota(data.begin(), data.end(), 0.0f);

  std::cout << "\n--- Use `process` to process each n samples at once.\n";
  std::array<float, 8> input;
  for (size_t i = 0; i < data.size(); i += 8) {
    for (size_t j = 0; j < input.size(); ++j) input[j] = data[i + j];
    std::cout << filter.process(input) << "\n";
  }

  std::cout << "\n--- Use `push` and `output` to get output for each sample.\n";
  filter.reset();
  for (size_t i = 0; i < data.size(); ++i) {
    filter.push(data[i]);
    if (i % 8 == 7) std::cout << filter.output() << "\n";
  }

  return 0;
}

{{
  {b0, b1, b2, a1, a2}, // Section 0
  {/*      ...     */}, // Section 1
  {/*      ...     */}, // Section 2
  // ...
}}