ピークホールドによるエンベロープ

import matplotlib.pyplot as plt
import numpy as np

def peakHoldForward(sig, holdtime, reset=0):
    """
    holdtime の単位はサンプル数。
    最後に最大値が更新されてから holdtime サンプル後に hold の値をリセットする。
    """
    out = np.empty(len(sig))
    hold = reset
    counter = 0
    for i in range(len(sig)):
        if counter > 0:  # 注意: 比較に >= を使うとホールド時間が 1 サンプル延びる。
            counter -= 1
        else:
            hold = reset
        if hold <= sig[i]:  # 注意: カウンタをリセットしたいので <= で比較。
            hold = sig[i]
            counter = holdtime
        out[i] = hold
    return out

rng = np.random.default_rng(0)

sig = rng.normal(0, 1, 256)
absed = np.abs(sig)
holdF = peakHoldForward(absed, 32, 0)

plt.plot(sig, label="Input")
plt.plot(absed, label="Absed")
plt.plot(holdF, label="Forward Hold")
plt.show()

def makeTriangleFir(delay):
    """delay は 2 より大きい整数。"""
    fir = np.interp(
        np.linspace(0, 1, delay + 1),
        [0, 0.5, 1],
        [0, 1, 0],
    )
    return fir / np.sum(fir)

from collections import deque

def idealPeakHoldNaive(sig, holdtime: int):
    out = np.empty(len(sig))
    buffer = deque([0 for _ in range(holdtime)])
    counter = 0
    for i in range(len(sig)):
        buffer.append(abs(sig[i]))
        buffer.popleft()
        out[i] = max(buffer)
    return out

from collections import deque

def idealPeakHoldFast(sig, holdTime):
    """
    sig: 入力信号。
    holdTime: ホールド時間。単位はサンプル数。
    """
    out = np.empty(len(sig))                     # 出力信号。
    buffer = deque([0 for _ in range(holdTime)]) # ディレイのバッファ。
    hold = deque([])                             # ホールド値のキュー。
    for i in range(len(sig)):
        x0 = sig[i]

        if len(hold) > 0: # 手順 1
            idx = len(hold) - 1
            while idx >= 0:
                if hold[idx] < x0:
                    hold.pop()
                else:
                    break
                idx -= 1

        hold.append(x0)   # 手順 2

        buffer.append(x0) # 手順 3
        delayOut = buffer.popleft()

        if len(hold) > 0 and delayOut == hold[0]: # 手順 4
            hold.popleft()

        out[i] = hold[0] if len(hold) > 0 else 0 # 手順 5

    return out

def nextTime(rng, rate):
    """
    ポアソン過程。
    もし単位時間当たり N 回のイベントが起こるなら `rate` を `1/N` に設定する。
    """
    return -np.log(1.0 - rng.uniform(0, 1)) / rate

def pulseNoise(rng, rate, length):
    out = np.zeros(length)
    time = 0
    while time < length:
        time += nextTime(rng, rate)

        t1 = int(time)
        t2 = t1 + 1
        amp = rng.uniform(0, 1)
        if t1 < length:
            out[t1] = amp * (t2 - time)
        if t2 < length:
            out[t2] = amp * (time - t1)
    return out

from multiprocessing import Pool

def verifyHoldEnvelope(sig, env, name):
    condition = (sig - env) > 8 * np.finfo(np.float64).eps
    if np.any(condition):
        idx = np.where(condition)
        diff = np.take(sig, idx) - np.take(env, idx)
        print(f"Test failed: {name}", diff, sep="\n")

def jobTestIdealHoldRandom(seed):
    hold = 32   # ホールド時間 (サンプル数)
    delay = 32  # フィルタの遅延 (サンプル数)
    triangleFir = makeTriangleFir(delay)

    rng = np.random.default_rng(seed)
    sig = pulseNoise(rng, 1 / 8, 48000)

    naive = idealPeakHoldNaive(sig, hold)
    fast = idealPeakHoldFast(sig, hold)

    smoothedNaive = signal.lfilter(triangleFir, 1, naive)
    smoothedFast = signal.lfilter(triangleFir, 1, fast)

    delayed = np.hstack((np.zeros(hold), sig[:-hold]))

    verifyHoldEnvelope(delayed, smoothedNaive, f"Naive, {seed}, {delay}")
    verifyHoldEnvelope(delayed, smoothedFast, f"Fast, {seed}, {delay}")

def testIdealHoldRandom(nTest=1024):
    with Pool() as pool:
        seeds = [i for i in range(nTest)]
        idx = 0
        for result in pool.imap_unordered(jobTestIdealHoldRandom, seeds):
            print(f"\r{idx}", end="")  # 処理が終わったシード値を表示。
            idx += 1

if __name__ == "__main__":
    testIdealHoldRandom()

#include <algorithm>
#include <iostream>
#include <random>
#include <vector>

// 整数サンプルのディレイ。
template<typename Sample> class IntDelay {
public:
  std::vector<Sample> buf;
  size_t wptr = 0;
  size_t rptr = 0;

  IntDelay(size_t size = 65536) : buf(size) {}

  void resize(size_t size)
  {
    buf.resize(size);
    wptr = 0;
    rptr = 0;
  }

  void reset() { std::fill(buf.begin(), buf.end(), Sample(0)); }

  void setFrames(size_t delayFrames)
  {
    if (delayFrames >= buf.size()) delayFrames = buf.size();
    rptr = wptr - delayFrames;
    if (rptr >= buf.size()) rptr += buf.size(); // Unsigned overflow case.
  }

  Sample process(Sample input)
  {
    if (++wptr >= buf.size()) wptr -= buf.size();
    buf[wptr] = input;

    if (++rptr >= buf.size()) rptr -= buf.size();
    return buf[rptr];
  }
};

// メモリの確保と解放を減らした std::deque の代用データ構造。
template<typename T> struct RingQueue {
  std::vector<T> buf;
  size_t wptr = 0;
  size_t rptr = 0;

  void resize(size_t size) { buf.resize(size); }

  void reset(T value = 0)
  {
    std::fill(buf.begin(), buf.end(), value);
    wptr = 0;
    rptr = 0;
  }

  inline size_t size()
  {
    auto sz = wptr - rptr;
    if (sz >= buf.size()) sz += buf.size(); // Unsigned overflow case.
    return sz;
  }

  inline bool empty() { return wptr == rptr; }

  T &front() { return buf[increment(rptr)]; }
  T &back() { return buf[wptr]; }

  inline size_t increment(size_t idx)
  {
    if (++idx >= buf.size()) idx -= buf.size();
    return idx;
  }

  inline size_t decrement(size_t idx)
  {
    if (--idx >= buf.size()) idx += buf.size(); // Unsigned overflow case.
    return idx;
  }

  void push_back(T value)
  {
    wptr = increment(wptr);
    buf[wptr] = value;
  }

  T pop_front()
  {
    rptr = increment(rptr);
    return buf[rptr];
  }

  T pop_back()
  {
    wptr = decrement(wptr);
    return buf[wptr];
  }
};

/*
理想的なホールドの実装。
- `setFrames(0)` のとき出力はすべて 0 。
- `setFrames(1)` のとき入力をバイパス。
*/
template<typename Sample> struct PeakHold {
  Sample neutral = 0;
  IntDelay<Sample> delay;
  RingQueue<Sample> hold;

  PeakHold(size_t size = 65536) {
    resize(size);
    setFrames(1);
  }

  void resize(size_t size)
  {
    delay.resize(size);
    hold.resize(size);
  }

  void reset()
  {
    delay.reset();
    hold.reset(neutral);
  }

  void setFrames(size_t frames) { delay.setFrames(frames); }

  Sample process(Sample x0)
  {
    while (!queue.empty()) {
      if (queue.back() >= x0) break;
      queue.pop_back();
    }
    queue.push_back(x0);
    if (delay.process(x0) == queue.front()) queue.pop_front();
    return queue.front();
  }
};

int main()
{
  // PeakHold の使用例。
  constexpr size_t length = 64;

  PeakHold<float> hold;
  hold.setFrames(4);

  std::minstd_rand rng(0);
  std::uniform_real_distribution<float> dist(0.0f, 1.0f);
  for (size_t idx = 0; idx < length; ++idx) {
    auto input = dist(rng);
    std::cout << input << ", " << hold.process(input) << "\n";
  }

  return 0;
}

def peakHoldBackward(sig, holdtime, reset=0):
    out = np.empty(len(sig))
    hold = reset
    counter = 0
    for i in reversed(range(len(sig))): # ここだけ変更。
        if counter > 0:
            counter -= 1
        else:
            hold = reset
        if hold <= sig[i]:
            hold = sig[i]
            counter = holdtime
        out[i] = hold
    return out