Highly efficient streaming sliding window (for peak hold)

[Fire Sledge - Ohm Force]

I was coding almost the same thing (max over a sliding window) and came with a binary tree algorithm. While searching for possibly more efficient methods, I found this thread. I think my method is quite similar to the OP’s algorithm, maybe slightly simpler because the window length remains fixed during a run (it only extends from 0 to N when inserting the first samples).

My code is a template and can take any binary operator with associative and commutative properties (min/max, addition, multiplication, string concatenation…) Here it is, if anyone is interested:

SlidingOp.h
SlidingOp.hpp

Typical use:

Code: Select all

#include <algorithm>
#include "mfx/dsp/iir/SlidingOp.h"

struct MaxFtor { float operator () (float a, float b) { return std::max (a, b); } };
…
mfx::dsp::fir::SlidingOp <float, MaxFtor> sliding_max;
sliding_max.set_length (123);
…
output = sliding_max.process_sample (input);

[Z1202]

FWIW a couple of years ago I designed a true O(1) sliding max algorithm. It was based on an amortized O(1) sliding max algorithm published on stackoverflow. I didn't store the link, but I think it's this one: https://stackoverflow.com/questions/480 ... all-consta

Here is the implementation. Some headers may be missing, depending on your compiler. Don't ask me for details, I don't remember them

Code: Select all

#include<array>
#include<assert.h>

/*

writepos  scan    scanpos  scanend   outpos  outrange   scanrange   inrange
   0       --        8        5        1       1..4       4..8 <       0  
   1       8 7       7        5        2       2..4       5..8        0..1
   2       7 6       6        5        3       3..4       5..8        0..2
   3       6 5       5        5        4        4         5..8        0..3
   4       --        4        0        5       5..8       0..3 <       4  
   5       4 3       3        0        6       6..8       0..3        4..5
   6       3 2       2        0        7       7..8       0..3        4..6
   7       2 1       1        0        8        8         0..3        4..7
   8       1 0       0        0        0       0..4       0..3        4..8

writepos  scan    scanpos  scanend   outpos  outrange   scanrange   inrange
   0       --        7        4        1       1..3       4..7 <       0
   1       7 6       6        4        2       2..3       4..7        0..1
   2       6 5       5        4        3        3         4..7        0..2
   3       5 4       4        4        4       4..7       4..7        0..3
   4       --        3        0        5       5..7       0..3 <       4
   5       3 2       2        0        6       6..7       0..3        4..5
   6       2 1       1        0        7        7         0..3        4..6
   7       1 0       0        0        0       0..3       0..3        4..7

*/

template<class T>
class MovingMax
{
public:
    enum { N = 8 };
private:
    enum {
        SCANSTARTLOW = (N-1)/2,
        SCANSTARTHIGH = N-1,
        SCANSTARTXOR = SCANSTARTLOW ^ SCANSTARTHIGH
    };
public:
    MovingMax()
    {
        m_scanmax = m_inmax = std::numeric_limits<T>::lowest();
        m_data.fill(m_inmax);
        m_writepos = 0;
        m_scanend = 0;
        m_scanpos = 0;
        m_scanstart = SCANSTARTLOW;
    }
    T operator<<(T x)
    {
        if( --m_scanpos >= m_scanend )
        {
            m_inmax = std::max(m_inmax,x);
            m_data[m_scanpos] = std::max( m_data[m_scanpos], m_data[m_scanpos+1] );
        }
        else
        {
            m_scanmax = m_inmax;
            m_inmax = x;
            m_scanend = m_scanend ^ ((N+1)/2);
            m_scanstart = m_scanstart ^ SCANSTARTXOR;
            m_scanpos = m_scanstart;
        }
 
        m_data[m_writepos] = x;
        if( ++m_writepos >= N )
            m_writepos = 0;
        T outmax = m_data[m_writepos];
        T movingmax = std::max( m_inmax, std::max(m_scanmax,outmax) );
        return movingmax;
    }
private:
    int m_writepos;
    int m_scanpos;
    int m_scanend;
    int m_scanstart;
    T m_inmax;
    T m_scanmax;
    std::array<T,N> m_data;
};

int main()
{
    using M = MovingMax<float>;
    M m;
    std::array<float,1000> a;
    int seed = 10;
    float y = 0;
    for( int i = 0 ; i < (int)a.size() ; i++ )
    {
        seed = seed*1103515245 + 12345;
        float x = seed/std::pow(2.f,24);
        a[i] = x;
    }
    for( int i = 0 ; i < (int)a.size() ; i++ )
    {
        float f = m << a[i];
        const int N=M::N;
        int i1 = std::max(i-N+1,0);
        float f1 = a[i1];
        int n = 1;
        for( i1++ ; i1 <= i ; i1++, n++ )
            f1 = std::max(f1,a[i1]);
        assert( f == f1 );
        assert( i < N || n == N );
    }
}

[chipnix]

Interesting thread!
Let me summarise, there are two approaches:

#1 binary tree approach with O(log n)
#2 two stacks approach with amortised O(1)

So which one is better for live-audio? A rough estimate for 64 sample callback-buffer, and a 8820-sample sized window (0,2 sec @ 44100)

(this is not measuring, this is only a rough estimate)

#1 log2(8820) per sample * 64 buffer-size * 4 (because of complexity)
≈ 1678 theoretical time units

#2. (almost) worst case, the algorithm must iterate once over the entire buffer size during the callback
≈ 8820 theoretical time units (worst case, but much less amortised)

It looks the binary tree approach is the better choice, especially for bigger window sizes for live audio.
So much for the theory, do any of you have practical experience? Has anyone measured both algorithms?



PS: @Z1202
not really understanding your algorithm, but I will figure it out. If its really real O(1), if this is working correct this would be the best option

[chipnix]

FWIW a couple of years ago I designed a true O(1) sliding max algorithm. It was based on an amortized O(1) sliding max algorithm published on stackoverflow. I didn't store the link, but I think it's this one: https://stackoverflow.com/questions/480 ... all-consta

Here is the implementation. Some headers may be missing, depending on your compiler. Don't ask me for details, I don't remember them

Code: Select all

#include<array>
#include<assert.h>

/*

writepos  scan    scanpos  scanend   outpos  outrange   scanrange   inrange
   0       --        8        5        1       1..4       4..8 <       0  
   1       8 7       7        5        2       2..4       5..8        0..1
   2       7 6       6        5        3       3..4       5..8        0..2
   3       6 5       5        5        4        4         5..8        0..3
   4       --        4        0        5       5..8       0..3 <       4  
   5       4 3       3        0        6       6..8       0..3        4..5
   6       3 2       2        0        7       7..8       0..3        4..6
   7       2 1       1        0        8        8         0..3        4..7
   8       1 0       0        0        0       0..4       0..3        4..8

writepos  scan    scanpos  scanend   outpos  outrange   scanrange   inrange
   0       --        7        4        1       1..3       4..7 <       0
   1       7 6       6        4        2       2..3       4..7        0..1
   2       6 5       5        4        3        3         4..7        0..2
   3       5 4       4        4        4       4..7       4..7        0..3
   4       --        3        0        5       5..7       0..3 <       4
   5       3 2       2        0        6       6..7       0..3        4..5
   6       2 1       1        0        7        7         0..3        4..6
   7       1 0       0        0        0       0..3       0..3        4..7

*/

template<class T>
class MovingMax
{
public:
    enum { N = 8 };
private:
    enum {
        SCANSTARTLOW = (N-1)/2,
        SCANSTARTHIGH = N-1,
        SCANSTARTXOR = SCANSTARTLOW ^ SCANSTARTHIGH
    };
public:
    MovingMax()
    {
        m_scanmax = m_inmax = std::numeric_limits<T>::lowest();
        m_data.fill(m_inmax);
        m_writepos = 0;
        m_scanend = 0;
        m_scanpos = 0;
        m_scanstart = SCANSTARTLOW;
    }
    T operator<<(T x)
    {
        if( --m_scanpos >= m_scanend )
        {
            m_inmax = std::max(m_inmax,x);
            m_data[m_scanpos] = std::max( m_data[m_scanpos], m_data[m_scanpos+1] );
        }
        else
        {
            m_scanmax = m_inmax;
            m_inmax = x;
            m_scanend = m_scanend ^ ((N+1)/2);
            m_scanstart = m_scanstart ^ SCANSTARTXOR;
            m_scanpos = m_scanstart;
        }
 
        m_data[m_writepos] = x;
        if( ++m_writepos >= N )
            m_writepos = 0;
        T outmax = m_data[m_writepos];
        T movingmax = std::max( m_inmax, std::max(m_scanmax,outmax) );
        return movingmax;
    }
private:
    int m_writepos;
    int m_scanpos;
    int m_scanend;
    int m_scanstart;
    T m_inmax;
    T m_scanmax;
    std::array<T,N> m_data;
};

int main()
{
    using M = MovingMax<float>;
    M m;
    std::array<float,1000> a;
    int seed = 10;
    float y = 0;
    for( int i = 0 ; i < (int)a.size() ; i++ )
    {
        seed = seed*1103515245 + 12345;
        float x = seed/std::pow(2.f,24);
        a[i] = x;
    }
    for( int i = 0 ; i < (int)a.size() ; i++ )
    {
        float f = m << a[i];
        const int N=M::N;
        int i1 = std::max(i-N+1,0);
        float f1 = a[i1];
        int n = 1;
        for( i1++ ; i1 <= i ; i1++, n++ )
            f1 = std::max(f1,a[i1]);
        assert( f == f1 );
        assert( i < N || n == N );
    }
}

I checked it, and it seems to work fantastic!

[avasopht]

I've got a feeling I might have just scrapped my code and used whatever was posted in here.

But kudos to Z1202 for their algorithm.

I do have a little DSP in the pipeline, but I honestly can't say if or when I'd run a benchmark to compare. I just recall the log2 algorithm being fast enough.

Sorry for the necropost, ... it's been a while