Filter.h

Go to the documentation of this file.

/***********************************************************************************************************************
*                                                                                                                      *
* libscopehal                                                                                                          *
*                                                                                                                      *
* Copyright (c) 2012-2026 Andrew D. Zonenberg and contributors                                                         *
* All rights reserved.                                                                                                 *
*                                                                                                                      *
* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the     *
* following conditions are met:                                                                                        *
*                                                                                                                      *
*    * Redistributions of source code must retain the above copyright notice, this list of conditions, and the         *
*      following disclaimer.                                                                                           *
*                                                                                                                      *
*    * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the       *
*      following disclaimer in the documentation and/or other materials provided with the distribution.                *
*                                                                                                                      *
*    * Neither the name of the author nor the names of any contributors may be used to endorse or promote products     *
*      derived from this software without specific prior written permission.                                           *
*                                                                                                                      *
* THIS SOFTWARE IS PROVIDED BY THE AUTHORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED   *
* TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL *
* THE AUTHORS BE HELD LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES        *
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR       *
* BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE       *
* POSSIBILITY OF SUCH DAMAGE.                                                                                          *
*                                                                                                                      *
***********************************************************************************************************************/
 
#ifndef Filter_h
#define Filter_h
 
#include "OscilloscopeChannel.h"
#include "FlowGraphNode.h"
#include "KahanSummation.h"
 
class QueueHandle;
 
class HistogramConstants
{
public:
    uint32_t    size;
    uint32_t    bins;
    float       nmin;
    float       nmax;
};
 
class WaveformCacheKey
{
public:
 
    WaveformCacheKey()
    : m_wfm(nullptr)
    , m_rev(0)
    {}
 
    WaveformCacheKey(WaveformBase* wfm)
    : m_wfm(wfm)
    , m_rev(wfm->m_revision)
    {}
 
    bool operator==(WaveformBase* wfm)
    { return (m_wfm == wfm) && (m_rev == wfm->m_revision); }
 
    bool operator==(WaveformCacheKey wfm)
    { return (m_wfm == wfm.m_wfm) && (m_rev == wfm.m_rev); }
 
    bool operator!=(WaveformBase* wfm)
    { return (m_wfm != wfm) || (m_rev != wfm->m_revision); }
 
    bool operator!=(WaveformCacheKey wfm)
    { return (m_wfm != wfm.m_wfm) || (m_rev != wfm.m_rev); }
 
    WaveformBase* m_wfm;
 
    uint64_t m_rev;
};
 
class Filter    : public OscilloscopeChannel
{
public:
 
    // Construction and enumeration
 
    enum Category
    {
        CAT_ANALYSIS,       //Signal integrity analysis
        CAT_BUS,            //Buses
        CAT_CLOCK,          //Clock stuff
        CAT_MATH,           //Basic math functions
        CAT_MEASUREMENT,    //Measurement functions
        CAT_MEMORY,         //Memory buses
        CAT_SERIAL,         //Serial communications
        CAT_MISC,           //anything not otherwise categorized
        CAT_POWER,          //Power analysis
        CAT_RF,             //Frequency domain analysis (FFT etc) and other RF stuff
        CAT_GENERATION,     //Waveform generation and synthesis
        CAT_EXPORT,         //Waveform export
        CAT_OPTICAL,        //Optics
 
        CAT_COUNT           //current max number of categories (will increase over time as more are added)
    };
 
    Filter(
        const std::string& color,
        Category cat,
        Unit xunit = Unit::UNIT_FS);
    virtual ~Filter();
 
    static const std::set<Filter*>& GetAllInstances()
    { return m_filters; }
 
    static size_t GetNumInstances()
    { return m_filters.size(); }
 
    void HideFromList()
    {
        m_filters.erase(this);
        m_instanceCount[GetProtocolDisplayName()] --;
    }
 
    virtual void ClearStreams() override;
    virtual size_t AddStream(
        Unit yunit,
        const std::string& name,
        Stream::StreamType stype,
        uint8_t flags = 0) override;
 
    void AddProtocolStream(const std::string& name)
    { AddStream(Unit(Unit::UNIT_COUNTS), name, Stream::STREAM_TYPE_PROTOCOL); }
 
    void AddDigitalStream(const std::string& name)
    { AddStream(Unit(Unit::UNIT_COUNTS), name, Stream::STREAM_TYPE_DIGITAL); }
 
    // Name generation
 
    virtual void SetDefaultName();
 
    void UseDefaultName(bool use)
    {
        m_usingDefault = use;
        if(use)
            SetDefaultName();
    }
 
    bool IsUsingDefaultName()
    { return m_usingDefault; }
 
    // Reference counting
 
    virtual void AddRef() override;
    virtual void Release() override;
 
    // Accessors
 
    Category GetCategory()
    { return m_category; }
 
    virtual std::string GetProtocolDisplayName() =0;
 
public:
    virtual void ClearSweeps();
 
    // Vertical scaling
 
    virtual void AutoscaleVertical(size_t stream);
 
    virtual float GetVoltageRange(size_t stream) override;
    virtual void SetVoltageRange(float range, size_t stream) override;
 
    virtual float GetOffset(size_t stream) override;
    virtual void SetOffset(float offset, size_t stream) override;
 
protected:
 
    std::vector<float> m_ranges;
 
    std::vector<float> m_offsets;
 
public:
    // Serialization
 
    virtual YAML::Node SerializeConfiguration(IDTable& table) override;
 
    virtual void LoadParameters(const YAML::Node& node, IDTable& table) override;
    virtual void LoadInputs(const YAML::Node& node, IDTable& table) override;
 
    virtual bool ShouldPersistWaveform() override;
 
protected:
 
    Category m_category;
 
    bool m_usingDefault;
 
    bool VerifyAllInputsOK(bool allowEmpty = false);
    bool VerifyInputOK(size_t i, bool allowEmpty = false);
    bool VerifyAllInputsOKAndUniformAnalog();
    bool VerifyAllInputsOKAndSparseAnalog();
    bool VerifyAllInputsOKAndSparseDigital();
    bool VerifyAllInputsOKAndSparseOrUniformDigital();
 
public:
    static int64_t GetNextEventTimestamp(SparseWaveformBase* wfm, size_t i, size_t len, int64_t timestamp);
    static int64_t GetNextEventTimestamp(UniformWaveformBase* wfm, size_t i, size_t len, int64_t timestamp);
 
    static int64_t GetNextEventTimestamp(
        SparseWaveformBase* swfm, UniformWaveformBase* uwfm, size_t i, size_t len, int64_t timestamp)
    {
        if(swfm)
            return GetNextEventTimestamp(swfm, i, len, timestamp);
        else
            return GetNextEventTimestamp(uwfm, i, len, timestamp);
    }
 
    static void AdvanceToTimestamp(SparseWaveformBase* wfm, size_t& i, size_t len, int64_t timestamp);
    static void AdvanceToTimestamp(UniformWaveformBase* wfm, size_t& i, size_t len, int64_t timestamp);
 
    static void AdvanceToTimestamp(
        SparseWaveformBase* swfm, UniformWaveformBase* uwfm, size_t& i, size_t len, int64_t timestamp)
    {
        if(swfm)
            AdvanceToTimestamp(swfm, i, len, timestamp);
        else
            AdvanceToTimestamp(uwfm, i, len, timestamp);
    }
 
    static int64_t GetNextEventTimestampScaled(SparseWaveformBase* wfm, size_t i, size_t len, int64_t timestamp);
    static int64_t GetNextEventTimestampScaled(UniformWaveformBase* wfm, size_t i, size_t len, int64_t timestamp);
    static void AdvanceToTimestampScaled(SparseWaveformBase* wfm, size_t& i, size_t len, int64_t timestamp);
    static void AdvanceToTimestampScaled(UniformWaveformBase* wfm, size_t& i, size_t len, int64_t timestamp);
 
    static void AdvanceToTimestampScaled(
        SparseWaveformBase* swfm, UniformWaveformBase* uwfm, size_t& i, size_t len, int64_t timestamp)
    {
        if(swfm)
            AdvanceToTimestampScaled(swfm, i, len, timestamp);
        else
            AdvanceToTimestampScaled(uwfm, i, len, timestamp);
    }
 
    static int64_t GetNextEventTimestampScaled(
        SparseWaveformBase* swfm, UniformWaveformBase* uwfm, size_t i, size_t len, int64_t timestamp)
    {
        if(swfm)
            return GetNextEventTimestampScaled(swfm, i, len, timestamp);
        else
            return GetNextEventTimestampScaled(uwfm, i, len, timestamp);
    }
 
protected:
    UniformAnalogWaveform* SetupEmptyUniformAnalogOutputWaveform(WaveformBase* din, size_t stream, bool clear=true);
    SparseAnalogWaveform* SetupEmptySparseAnalogOutputWaveform(WaveformBase* din, size_t stream, bool clear=true);
    UniformDigitalWaveform* SetupEmptyUniformDigitalOutputWaveform(WaveformBase* din, size_t stream);
    SparseDigitalWaveform* SetupEmptySparseDigitalOutputWaveform(WaveformBase* din, size_t stream);
    SparseAnalogWaveform* SetupSparseOutputWaveform(SparseWaveformBase* din, size_t stream, size_t skipstart, size_t skipend);
    SparseDigitalWaveform* SetupSparseDigitalOutputWaveform(SparseWaveformBase* din, size_t stream, size_t skipstart, size_t skipend);
 
    template<class T>
    T* SetupEmptyWaveform(WaveformBase* din, size_t stream, bool clear = true)
    {
        //Create the waveform, but only if necessary
        auto cap = dynamic_cast<T*>(GetData(stream));
        if(cap == nullptr)
        {
            cap = new T;
            SetData(cap, stream);
        }
 
        //Copy configuration
        if(din != nullptr)
        {
            cap->m_startTimestamp       = din->m_startTimestamp;
            cap->m_startFemtoseconds    = din->m_startFemtoseconds;
            cap->m_triggerPhase         = din->m_triggerPhase;
            cap->m_timescale            = din->m_timescale;
        }
 
        //Bump rev number
        cap->m_revision ++;
 
        //Clear output
        if(clear)
            cap->clear();
 
        return cap;
    }
 
public:
    //Helpers for sub-sample interpolation
 
    template<class T>
    __attribute__((noinline))
    static float InterpolateTime(T* cap, size_t a, float voltage)
    {
        AssertTypeIsAnalogWaveform(cap);
 
        //If the voltage isn't between the two points, abort
        float fa = cap->m_samples[a];
        float fb = cap->m_samples[a+1];
        bool ag = (fa > voltage);
        bool bg = (fb > voltage);
        if( (ag && bg) || (!ag && !bg) )
            return 0;
 
        //no need to divide by time, sample spacing is normalized to 1 timebase unit
        float slope = (fb - fa);
        float delta = voltage - fa;
        return delta / slope;
    }
 
    static float InterpolateTime(SparseAnalogWaveform* s, UniformAnalogWaveform* u, size_t a, float voltage)
    {
        if(s)
            return InterpolateTime(s, a, voltage);
        else
            return InterpolateTime(u, a, voltage);
    }
 
    static float InterpolateTime(UniformAnalogWaveform* p, UniformAnalogWaveform* n, size_t a, float voltage);
    static float InterpolateTime(SparseAnalogWaveform* p, SparseAnalogWaveform* n, size_t a, float voltage);
 
    static float InterpolateTime(
        SparseAnalogWaveform* sp,
        UniformAnalogWaveform* up,
        SparseAnalogWaveform* sn,
        UniformAnalogWaveform* un,
        size_t a, float voltage)
    {
        if(sp)
            return InterpolateTime(sp, sn, a, voltage);
        else
            return InterpolateTime(up, un, a, voltage);
    }
 
    static float InterpolateValue(SparseAnalogWaveform* cap, size_t index, float frac_ticks);
    static float InterpolateValue(UniformAnalogWaveform* cap, size_t index, float frac_ticks);
 
    //Helpers for more complex measurements
    //TODO: create some process for caching this so we don't waste CPU time
 
    template<class T>
    __attribute__((noinline))
    static void GetMinMaxVoltage(T* cap, float& vmin, float& vmax)
    {
        AssertTypeIsAnalogWaveform(cap);
 
        vmin = FLT_MAX;
        vmax = FLT_MIN;
        for(float f : cap->m_samples)
        {
            if(f < vmin)
                vmin = f;
            if(f > vmax)
                vmax = f;
        }
    }
 
    template<class T>
    __attribute__((noinline))
    static void GetMinMaxVoltage(
        vk::raii::CommandBuffer& cmdBuf,
        std::shared_ptr<QueueHandle> queue,
        ComputePipeline& minmaxPipeline,
        AcceleratorBuffer<float>& scratchMin,
        AcceleratorBuffer<float>& scratchMax,
        T* cap,
        float& vmin,
        float& vmax
        )
    {
        AssertTypeIsAnalogWaveform(cap);
 
        //GPU side min/max
        const uint32_t nthreads = 4096;
        const uint32_t threadsPerBlock = 64;
        scratchMin.resize(nthreads);
        scratchMax.resize(nthreads);
 
        cmdBuf.begin({});
        minmaxPipeline.BindBufferNonblocking(0, cap->m_samples, cmdBuf);
        minmaxPipeline.BindBufferNonblocking(1, scratchMin, cmdBuf);
        minmaxPipeline.BindBufferNonblocking(2, scratchMax, cmdBuf);
        minmaxPipeline.Dispatch(cmdBuf, (uint32_t)cap->size(), nthreads / threadsPerBlock);
 
        scratchMin.MarkModifiedFromGpu();
        scratchMax.MarkModifiedFromGpu();
 
        scratchMin.PrepareForCpuAccessNonblocking(cmdBuf);
        scratchMax.PrepareForCpuAccessNonblocking(cmdBuf);
 
        cmdBuf.end();
        queue->SubmitAndBlock(cmdBuf);
 
        //Final reduction on CPU (TODO: faster to do second shader invocation and read one value??)
        vmin = scratchMin[0];
        vmax = scratchMax[0];
        for(uint32_t i=1; i<nthreads; i++)
        {
            vmin = std::min(vmin, scratchMin[i]);
            vmax = std::max(vmax, scratchMax[i]);
        }
    }
 
    template<class T>
    __attribute__((noinline))
    static float GetMinVoltage(T* cap)
    {
        AssertTypeIsAnalogWaveform(cap);
 
        //Loop over samples and find the minimum
        float tmp = FLT_MAX;
        for(float f : cap->m_samples)
        {
            if(f < tmp)
                tmp = f;
        }
        return tmp;
    }
 
    static float GetMinVoltage(SparseAnalogWaveform* s, UniformAnalogWaveform* u)
    {
        if(s)
            return GetMinVoltage(s);
        else
            return GetMinVoltage(u);
    }
 
    template<class T>
    __attribute__((noinline))
    static float GetMaxVoltage(T* cap)
    {
        AssertTypeIsAnalogWaveform(cap);
 
        //Loop over samples and find the maximum
        float tmp = -FLT_MAX;
        for(float f : cap->m_samples)
        {
            if(f > tmp)
                tmp = f;
        }
        return tmp;
    }
 
    static float GetMaxVoltage(SparseAnalogWaveform* s, UniformAnalogWaveform* u)
    {
        if(s)
            return GetMaxVoltage(s);
        else
            return GetMaxVoltage(u);
    }
 
    static void GetBaseAndTopVoltage(
        vk::raii::CommandBuffer& cmdBuf,
        std::shared_ptr<QueueHandle> queue,
        ComputePipeline& minmaxPipeline,
        std::shared_ptr<ComputePipeline>& histogramPipeline,
        AcceleratorBuffer<float>& scratchMin,
        AcceleratorBuffer<float>& scratchMax,
        AcceleratorBuffer<uint64_t>& hist,
        SparseAnalogWaveform* swfm,
        UniformAnalogWaveform* uwfm,
        float& base,
        float& top)
    {
        if(swfm)
        {
            GetBaseAndTopVoltage(
                cmdBuf,
                queue,
                minmaxPipeline,
                histogramPipeline,
                scratchMin,
                scratchMax,
                hist,
                swfm,
                base,
                top);
        }
        else if(uwfm)
        {
            GetBaseAndTopVoltage(
                cmdBuf,
                queue,
                minmaxPipeline,
                histogramPipeline,
                scratchMin,
                scratchMax,
                hist,
                uwfm,
                base,
                top);
        }
        else
        {
            LogError("Filter::GetBaseAndTopVoltage called with both waveforms null\n");
            base = 0;
            top = 0;
        }
    }
 
    template<class T>
    __attribute__((noinline))
    static void GetBaseAndTopVoltage(
        vk::raii::CommandBuffer& cmdBuf,
        std::shared_ptr<QueueHandle> queue,
        ComputePipeline& minmaxPipeline,
        std::shared_ptr<ComputePipeline>& histogramPipeline,
        AcceleratorBuffer<float>& scratchMin,
        AcceleratorBuffer<float>& scratchMax,
        AcceleratorBuffer<uint64_t>& hist,
 
        T* cap, float& base, float& top)
    {
        AssertTypeIsAnalogWaveform(cap);
 
        //GPU side min/max calculation
        float vmin;
        float vmax;
        GetMinMaxVoltage(cmdBuf, queue, minmaxPipeline, scratchMin, scratchMax, cap, vmin, vmax);
 
        //Do the histogram calculation on the GPU if possible
        float delta = vmax - vmin;
        const int nbins = 100;
        if(g_hasShaderInt64 && g_hasShaderAtomicInt64 && histogramPipeline)
        {
            MakeHistogram(cmdBuf, queue, *histogramPipeline, cap, hist, vmin, vmax, nbins);
            hist.PrepareForCpuAccess();
        }
        else
        {
            //CPU fallback
            hist.resize(nbins);
            hist.PrepareForCpuAccess();
 
            auto tmp = MakeHistogram(cap, vmin, vmax, nbins);
            for(int i=0; i<nbins; i++)
                hist[i] = tmp[i];
 
            hist.MarkModifiedFromCpu();
        }
 
        //Find the highest peak in the first quarter of the histogram
        size_t binval = 0;
        int idx = 0;
        for(int i=0; i<(nbins/4); i++)
        {
            if(hist[i] > binval)
            {
                binval = hist[i];
                idx = i;
            }
        }
 
        float fbin = (idx + 0.5f)/nbins;
        base = fbin*delta + vmin;
 
        //Find the highest peak in the third quarter of the histogram
        binval = 0;
        idx = 0;
        for(int i=(nbins*3)/4; i<nbins; i++)
        {
            if(hist[i] > binval)
            {
                binval = hist[i];
                idx = i;
            }
        }
 
        fbin = (idx + 0.5f)/nbins;
        top = fbin*delta + vmin;
    }
 
    template<class T>
    __attribute__((noinline))
    static float GetBaseVoltage(T* cap)
    {
        AssertTypeIsAnalogWaveform(cap);
 
        float vmin;
        float vmax;
        GetMinMaxVoltage(cap, vmin, vmax);
 
        float delta = vmax - vmin;
        const int nbins = 100;
        auto hist = MakeHistogram(cap, vmin, vmax, nbins);
 
        //Find the highest peak in the first quarter of the histogram
        size_t binval = 0;
        int idx = 0;
        for(int i=0; i<(nbins/4); i++)
        {
            if(hist[i] > binval)
            {
                binval = hist[i];
                idx = i;
            }
        }
 
        float fbin = (idx + 0.5f)/nbins;
        return fbin*delta + vmin;
    }
 
    static float GetBaseVoltage(SparseAnalogWaveform* swfm, UniformAnalogWaveform* uwfm)
    {
        if(swfm)
            return GetBaseVoltage(swfm);
        else
            return GetBaseVoltage(uwfm);
    }
 
    template<class T>
    __attribute__((noinline))
    static float GetTopVoltage(T* cap)
    {
        AssertTypeIsAnalogWaveform(cap);
 
        float vmin;
        float vmax;
        GetMinMaxVoltage(cap, vmin, vmax);
 
        float delta = vmax - vmin;
        const int nbins = 100;
        auto hist = MakeHistogram(cap, vmin, vmax, nbins);
 
        //Find the highest peak in the third quarter of the histogram
        size_t binval = 0;
        int idx = 0;
        for(int i=(nbins*3)/4; i<nbins; i++)
        {
            if(hist[i] > binval)
            {
                binval = hist[i];
                idx = i;
            }
        }
 
        float fbin = (idx + 0.5f)/nbins;
        return fbin*delta + vmin;
    }
 
    static float GetTopVoltage(SparseAnalogWaveform* swfm, UniformAnalogWaveform* uwfm)
    {
        if(swfm)
            return GetTopVoltage(swfm);
        else
            return GetTopVoltage(uwfm);
    }
 
    template<class T>
    __attribute__((noinline))
    static float GetAvgVoltage(T* cap)
    {
        AssertTypeIsAnalogWaveform(cap);
 
        //Loop over samples and find the average
        KahanSummation ksum;
        for(float f : cap->m_samples)
            ksum += f;
        return ksum.GetSum() / cap->m_samples.size();
    }
 
    static float GetAvgVoltage(SparseAnalogWaveform* swfm, UniformAnalogWaveform* uwfm)
    {
        if(swfm)
            return GetAvgVoltage(swfm);
        else
            return GetAvgVoltage(uwfm);
    }
 
    template<class T>
    __attribute__((noinline))
    static std::vector<size_t> MakeHistogram(T* cap, float low, float high, size_t bins)
    {
        AssertTypeIsAnalogWaveform(cap);
 
        std::vector<size_t> ret;
        for(size_t i=0; i<bins; i++)
            ret.push_back(0);
 
        //Early out if we have zero span
        if(bins == 0)
            return ret;
 
        float delta = high-low;
 
        for(float v : cap->m_samples)
        {
            float fbin = (v-low) / delta;
            size_t bin = floor(fbin * bins);
            if(fbin < 0)
                bin = 0;
            else
                bin = std::min(bin, bins-1);
            ret[bin] ++;
        }
 
        return ret;
    }
 
    static std::vector<size_t> MakeHistogram(
        SparseAnalogWaveform* s, UniformAnalogWaveform* u, float low, float high, size_t bins)
    {
        if(s)
            return MakeHistogram(s, low, high, bins);
        else
            return MakeHistogram(u, low, high, bins);
    }
 
    template<class T>
    __attribute__((noinline))
    static void MakeHistogram(
        vk::raii::CommandBuffer& cmdBuf,
        std::shared_ptr<QueueHandle> queue,
        ComputePipeline& histogramPipeline,
        T* cap,
        AcceleratorBuffer<uint64_t>& hist,
        float low,
        float high,
        size_t bins)
    {
        AssertTypeIsAnalogWaveform(cap);
 
        //Early out if we have zero span
        if(bins == 0)
            return;
 
        const uint32_t nthreads = 4096;
        const uint32_t threadsPerBlock = 8;
 
        //Fill the input histogram with zeroes CPU side for now
        //TODO: GPU side fill
        hist.resize(bins);
        hist.PrepareForCpuAccess();
        memset(hist.GetCpuPointer(), 0, bins*sizeof(int64_t));
        hist.MarkModifiedFromCpu();
 
        //Push constants
        HistogramConstants cfg;
        cfg.size = cap->size();
        cfg.bins = bins;
        cfg.nmin = low;
        cfg.nmax = high;
 
        cmdBuf.begin({});
        histogramPipeline.BindBufferNonblocking(0, cap->m_samples, cmdBuf);
        histogramPipeline.BindBufferNonblocking(1, hist, cmdBuf);
        histogramPipeline.Dispatch(cmdBuf, cfg, nthreads / threadsPerBlock);
 
        hist.MarkModifiedFromGpu();
 
        cmdBuf.end();
        queue->SubmitAndBlock(cmdBuf);
    }
 
    template<class T>
    __attribute__((noinline))
    static std::vector<size_t> MakeHistogramClipped(T* cap, float low, float high, size_t bins)
    {
        AssertTypeIsAnalogWaveform(cap);
 
        std::vector<size_t> ret;
        for(size_t i=0; i<bins; i++)
            ret.push_back(0);
 
        //Early out if we have zero span
        if(bins == 0)
            return ret;
 
        float delta = high-low;
 
        for(float v : cap->m_samples)
        {
            float fbin = (v-low) / delta;
            // must cast through a signed int type to avoid UB (e.g. saturates to 0 on arm64) [conv.fpint]
            size_t bin = static_cast<ssize_t>(floor(fbin * bins));
            if(bin >= bins) //negative values wrap to huge positive and get caught here
                continue;
            ret[bin] ++;
        }
 
        return ret;
    }
 
    template<class T, class R, class S>
    __attribute__((noinline))
    static void SampleOnAnyEdges(T* data, R* clock, SparseWaveform<S>& samples, bool cpuOnly = true)
    {
        //Compile-time check to make sure inputs are correct types
        AssertTypeIsDigitalWaveform(clock);
        AssertSampleTypesAreSame(data, &samples);
 
        samples.clear();
        if(cpuOnly)
            samples.SetGpuAccessHint(AcceleratorBuffer<S>::HINT_NEVER); //assume we're being used as part of a CPU-side filter
        clock->PrepareForCpuAccess();
        data->PrepareForCpuAccess();
        samples.PrepareForCpuAccess();
 
        //TODO: split up into blocks and multithread?
        //TODO: AVX vcompress?
        size_t len = clock->size();
        size_t dlen = data->size();
 
        //If the clock is sparse, assume it probably has edges on every sample and allocate that much buffer to start
        //(we might overallocate here but it'll be a lot faster)
        if(dynamic_cast<SparseDigitalWaveform*>(clock) != nullptr)
        {
            //Allocate exactly enough space
            samples.Resize(clock->size());
 
            size_t ndata = 0;
            size_t nout = 0;
            for(size_t i=1; i<len; i++)
            {
                //Throw away clock samples until we find an edge
                if(clock->m_samples[i] == clock->m_samples[i-1])
                    continue;
 
                //Throw away data samples until the data is synced with us
                int64_t clkstart = GetOffsetScaled(clock, i);
                while( (ndata+1 < dlen) && (GetOffsetScaled(data, ndata+1) < clkstart) )
                    ndata ++;
                if(ndata >= dlen)
                    break;
 
                //Add the new sample
                samples.m_offsets[nout] = clkstart;
                samples.m_samples[nout] = data->m_samples[ndata];
                nout ++;
            }
            samples.Resize(nout);
        }
        else
        {
            samples.Reserve(1 * 1024 * 1024);   //preallocate 1 MB sample buffer to avoid lots of reallocation when small
                                                //if it's smaller than this, we won't waste a lot of memory
 
            size_t ndata = 0;
            for(size_t i=1; i<len; i++)
            {
                //Throw away clock samples until we find an edge
                if(clock->m_samples[i] == clock->m_samples[i-1])
                    continue;
 
                //Throw away data samples until the data is synced with us
                int64_t clkstart = GetOffsetScaled(clock, i);
                while( (ndata+1 < dlen) && (GetOffsetScaled(data, ndata+1) < clkstart) )
                    ndata ++;
                if(ndata >= dlen)
                    break;
 
                //Add the new sample
                samples.m_offsets.push_back(clkstart);
                samples.m_samples.push_back(data->m_samples[ndata]);
            }
        }
        samples.MarkModifiedFromCpu();
 
        //Compute sample durations
        #ifdef __x86_64__
        if(g_hasAvx2)
            FillDurationsAVX2(samples);
        else
        #endif
            FillDurationsGeneric(samples);
 
        samples.MarkModifiedFromCpu();
    }
 
    template<class T>
    __attribute__((noinline))
    static void SampleOnAnyEdgesBase(WaveformBase* data, WaveformBase* clock, SparseWaveform<T>& samples)
    {
        data->PrepareForCpuAccess();
        clock->PrepareForCpuAccess();
        samples.PrepareForCpuAccess();
 
        auto udata = dynamic_cast<UniformWaveform<T>*>(data);
        auto sdata = dynamic_cast<SparseWaveform<T>*>(data);
 
        auto uclock = dynamic_cast<UniformDigitalWaveform*>(clock);
        auto sclock = dynamic_cast<SparseDigitalWaveform*>(clock);
 
        if(udata && uclock)
            SampleOnAnyEdges(udata, uclock, samples);
        else if(udata && sclock)
            SampleOnAnyEdges(udata, sclock, samples);
        else if(sdata && sclock)
            SampleOnAnyEdges(sdata, sclock, samples);
        else if(sdata && uclock)
            SampleOnAnyEdges(sdata, uclock, samples);
    }
 
    template<class T, class R, class S>
    __attribute__((noinline))
    static void SampleOnRisingEdges(T* data, R* clock, SparseWaveform<S>& samples)
    {
        //Compile-time check to make sure inputs are correct types
        AssertTypeIsDigitalWaveform(clock);
        AssertTypeIsSparseWaveform(&samples);
        AssertSampleTypesAreSame(data, &samples);
 
        samples.clear();
        samples.SetGpuAccessHint(AcceleratorBuffer<S>::HINT_NEVER); //assume we're being used as part of a CPU-side filter
        samples.Reserve(1 * 1024 * 1024);   //preallocate 1 MB sample buffer to avoid lots of reallocation when small
                                            //if it's smaller than this, we won't waste a lot of memory
 
        //TODO: split up into blocks and multithread?
        //TODO: AVX vcompress?
 
        size_t len = clock->size();
        size_t dlen = data->size();
 
        size_t ndata = 0;
        for(size_t i=1; i<len; i++)
        {
            //Throw away clock samples until we find a rising edge
            if(!(clock->m_samples[i] && !clock->m_samples[i-1]))
                continue;
 
            //Throw away data samples until the data is synced with us
            int64_t clkstart = GetOffsetScaled(clock, i);
            while( (ndata+1 < dlen) && (GetOffsetScaled(data, ndata+1) < clkstart) )
                ndata ++;
            if(ndata >= dlen)
                break;
 
            //Add the new sample
            samples.m_offsets.push_back(clkstart);
            samples.m_samples.push_back(data->m_samples[ndata]);
        }
 
        //Compute sample durations
        #ifdef __x86_64__
        if(g_hasAvx2)
            FillDurationsAVX2(samples);
        else
        #endif
            FillDurationsGeneric(samples);
 
        samples.MarkModifiedFromCpu();
    }
 
    template<class T>
    __attribute__((noinline))
    static void SampleOnRisingEdgesBase(WaveformBase* data, WaveformBase* clock, SparseWaveform<T>& samples)
    {
        data->PrepareForCpuAccess();
        clock->PrepareForCpuAccess();
        samples.PrepareForCpuAccess();
 
        auto udata = dynamic_cast<UniformWaveform<T>*>(data);
        auto sdata = dynamic_cast<SparseWaveform<T>*>(data);
 
        auto uclock = dynamic_cast<UniformDigitalWaveform*>(clock);
        auto sclock = dynamic_cast<SparseDigitalWaveform*>(clock);
 
        if(udata && uclock)
            SampleOnRisingEdges(udata, uclock, samples);
        else if(udata && sclock)
            SampleOnRisingEdges(udata, sclock, samples);
        else if(sdata && sclock)
            SampleOnRisingEdges(sdata, sclock, samples);
        else if(sdata && uclock)
            SampleOnRisingEdges(sdata, uclock, samples);
    }
 
    template<class T, class R, class S>
    __attribute__((noinline))
    static void SampleOnFallingEdges(T* data, R* clock, SparseWaveform<S>& samples)
    {
        //Compile-time check to make sure inputs are correct types
        AssertTypeIsDigitalWaveform(clock);
        AssertTypeIsSparseWaveform(&samples);
        AssertSampleTypesAreSame(data, &samples);
 
        samples.clear();
        samples.SetGpuAccessHint(AcceleratorBuffer<S>::HINT_NEVER); //assume we're being used as part of a CPU-side filter
        samples.Reserve(1 * 1024 * 1024);   //preallocate 1 MB sample buffer to avoid lots of reallocation when small
                                            //if it's smaller than this, we won't waste a lot of memory
 
        //TODO: split up into blocks and multithread?
        //TODO: AVX vcompress?
 
        size_t len = clock->size();
        size_t dlen = data->size();
 
        size_t ndata = 0;
        for(size_t i=1; i<len; i++)
        {
            //Throw away clock samples until we find a falling edge
            if(!(!clock->m_samples[i] && clock->m_samples[i-1]))
                continue;
 
            //Throw away data samples until the data is synced with us
            int64_t clkstart = GetOffsetScaled(clock, i);
            while( (ndata+1 < dlen) && (GetOffsetScaled(data, ndata+1) < clkstart) )
                ndata ++;
            if(ndata >= dlen)
                break;
 
            //Add the new sample
            samples.m_offsets.push_back(clkstart);
            samples.m_samples.push_back(data->m_samples[ndata]);
        }
 
        //Compute sample durations
        #ifdef __x86_64__
        if(g_hasAvx2)
            FillDurationsAVX2(samples);
        else
        #endif
            FillDurationsGeneric(samples);
 
        samples.MarkModifiedFromCpu();
    }
 
    template<class T, class R>
    __attribute__((noinline))
    static void SampleOnAnyEdgesWithInterpolation(T* data, R* clock, SparseAnalogWaveform& samples)
    {
        //Compile-time check to make sure inputs are correct types
        AssertTypeIsAnalogWaveform(data);
        AssertTypeIsDigitalWaveform(clock);
 
        samples.clear();
        samples.SetGpuAccessHint(AcceleratorBuffer<float>::HINT_NEVER); //assume we're being used as part of a CPU-side filter
        samples.Reserve(1 * 1024 * 1024);   //preallocate 1 MB sample buffer to avoid lots of reallocation when small
                                            //if it's smaller than this, we won't waste a lot of memory
 
        //TODO: split up into blocks and multithread?
        //TODO: AVX vcompress
 
        size_t len = clock->size();
        size_t dlen = data->size();
 
        size_t ndata = 0;
        for(size_t i=1; i<len; i++)
        {
            //Throw away clock samples until we find an edge
            if(clock->m_samples[i] == clock->m_samples[i-1])
                continue;
 
            //Throw away data samples until the data is synced with us
            int64_t clkstart = GetOffsetScaled(clock, i);
            while( (ndata+1 < dlen) && (GetOffsetScaled(data, ndata+1) < clkstart) )
                ndata ++;
            if(ndata >= dlen)
                break;
 
            //Find the fractional position of the clock edge
            int64_t tsample = GetOffsetScaled(data, ndata);
            int64_t delta = clkstart - tsample;
            float frac = delta * 1.0 / data->m_timescale;
 
            //Add the new sample
            samples.m_offsets.push_back(clkstart);
            samples.m_samples.push_back(InterpolateValue(data, ndata, frac));
        }
 
        //Compute sample durations
        #ifdef __x86_64__
        if(g_hasAvx2)
            FillDurationsAVX2(samples);
        else
        #endif
            FillDurationsGeneric(samples);
 
        samples.MarkModifiedFromCpu();
    }
 
    template<class T>
    __attribute__((noinline))
    static void SampleOnAnyEdgesBaseWithInterpolation(WaveformBase* data, WaveformBase* clock, SparseWaveform<T>& samples)
    {
        data->PrepareForCpuAccess();
        clock->PrepareForCpuAccess();
        samples.PrepareForCpuAccess();
 
        auto udata = dynamic_cast<UniformWaveform<T>*>(data);
        auto sdata = dynamic_cast<SparseWaveform<T>*>(data);
 
        auto uclock = dynamic_cast<UniformDigitalWaveform*>(clock);
        auto sclock = dynamic_cast<SparseDigitalWaveform*>(clock);
 
        if(udata && uclock)
            SampleOnAnyEdgesWithInterpolation(udata, uclock, samples);
        else if(udata && sclock)
            SampleOnAnyEdgesWithInterpolation(udata, sclock, samples);
        else if(sdata && sclock)
            SampleOnAnyEdgesWithInterpolation(sdata, sclock, samples);
        else if(sdata && uclock)
            SampleOnAnyEdgesWithInterpolation(sdata, uclock, samples);
    }
 
    template<class T>
    static void PrepareForCpuAccess(SparseWaveform<T>* s, UniformWaveform<T>* u)
    {
        if(s)
            s->PrepareForCpuAccess();
        else
            u->PrepareForCpuAccess();
    }
 
    template<class T>
    static void PrepareForGpuAccess(SparseWaveform<T>* s, UniformWaveform<T>* u)
    {
        if(s)
            s->PrepareForGpuAccess();
        else
            u->PrepareForGpuAccess();
    }
 
    static void FindRisingEdges(UniformAnalogWaveform* data, float threshold, std::vector<int64_t>& edges);
    static void FindRisingEdges(SparseAnalogWaveform* data, float threshold, std::vector<int64_t>& edges);
    static void FindZeroCrossings(SparseAnalogWaveform* data, float threshold, std::vector<int64_t>& edges);
    static void FindZeroCrossings(UniformAnalogWaveform* data, float threshold, std::vector<int64_t>& edges);
    static void FindZeroCrossings(UniformDigitalWaveform* data, std::vector<int64_t>& edges);
    static void FindZeroCrossings(SparseDigitalWaveform* data, std::vector<int64_t>& edges);
    static void FindRisingEdges(UniformDigitalWaveform* data, std::vector<int64_t>& edges);
    static void FindRisingEdges(SparseDigitalWaveform* data, std::vector<int64_t>& edges);
    static void FindFallingEdges(UniformDigitalWaveform* data, std::vector<int64_t>& edges);
    static void FindFallingEdges(SparseDigitalWaveform* data, std::vector<int64_t>& edges);
    static void FindPeaks(UniformAnalogWaveform* data, float peak_threshold, std::vector<int64_t>& peak_indices);
    static void FindPeaks(SparseAnalogWaveform* data, float peak_threshold, std::vector<int64_t>& peak_indices);
 
    static void FindZeroCrossingsBase(WaveformBase* data, float threshold, std::vector<int64_t>& edges)
    {
        auto udata = dynamic_cast<UniformAnalogWaveform*>(data);
        auto sdata = dynamic_cast<SparseAnalogWaveform*>(data);
 
        if(udata)
            FindZeroCrossings(udata, threshold, edges);
        else
            FindZeroCrossings(sdata, threshold, edges);
    }
 
    static void FindRisingEdges(
        SparseDigitalWaveform* sdata, UniformDigitalWaveform* udata, std::vector<int64_t>& edges)
    {
        if(sdata)
            FindRisingEdges(sdata, edges);
        else
            FindRisingEdges(udata, edges);
    }
 
    static void FindFallingEdges(
        SparseDigitalWaveform* sdata, UniformDigitalWaveform* udata, std::vector<int64_t>& edges)
    {
        if(sdata)
            FindFallingEdges(sdata, edges);
        else
            FindFallingEdges(udata, edges);
    }
 
    static void FindPeaks(
        SparseAnalogWaveform* sdata, UniformAnalogWaveform* udata, float peak_threshold, std::vector<int64_t>& peak_indices)
    {
        if(sdata)
            FindPeaks(sdata, peak_threshold, peak_indices);
        else
            FindPeaks(udata, peak_threshold, peak_indices);
    }
 
    static void FindZeroCrossings(
        SparseAnalogWaveform* sdata, UniformAnalogWaveform* udata, float threshold, std::vector<int64_t>& edges)
    {
        if(sdata)
            FindZeroCrossings(sdata, threshold, edges);
        else
            FindZeroCrossings(udata, threshold, edges);
    }
 
    static void FindZeroCrossings(
        SparseDigitalWaveform* sdata, UniformDigitalWaveform* udata, std::vector<int64_t>& edges)
    {
        if(sdata)
            FindZeroCrossings(sdata, edges);
        else
            FindZeroCrossings(udata, edges);
    }
 
    static void ClearAnalysisCache();
 
    enum FIRFilterType
    {
        FILTER_TYPE_LOWPASS,
        FILTER_TYPE_HIGHPASS,
        FILTER_TYPE_BANDPASS,
        FILTER_TYPE_NOTCH
    };
 
    static void CalculateFIRCoefficients(
        float fa,
        float fb,
        float stopbandAtten,
        FIRFilterType type,
        AcceleratorBuffer<float>& coefficients);
    static float Bessel(float x);
 
protected:
    //Helpers for sparse waveforms
    static void FillDurationsGeneric(SparseWaveformBase& wfm);
#ifdef __x86_64__
    static void FillDurationsAVX2(SparseWaveformBase& wfm);
#endif
 
public:
    sigc::signal<void()> signal_outputsChanged()
    { return m_outputsChangedSignal; }
 
protected:
    sigc::signal<void()> m_outputsChangedSignal;
 
    unsigned int m_instanceNum;
 
public:
    typedef Filter* (*CreateProcType)(const std::string&);
    static void DoAddDecoderClass(const std::string& name, CreateProcType proc);
 
    static void EnumProtocols(std::vector<std::string>& names);
    static Filter* CreateFilter(const std::string& protocol, const std::string& color = "#ffffff");
 
protected:
    //Class enumeration
    typedef std::map< std::string, CreateProcType > CreateMapType;
    static CreateMapType m_createprocs;
 
    //Object enumeration
    static std::set<Filter*> m_filters;
 
    //Instance naming
    static std::map<std::string, unsigned int> m_instanceCount;
 
    //Caching
    static std::mutex m_cacheMutex;
    static std::map<std::pair<WaveformBase*, float>, std::vector<int64_t> > m_zeroCrossingCache;
};
 
#define PROTOCOL_DECODER_INITPROC(T) \
    static Filter* CreateInstance(const std::string& color) \
    { \
        return new T(color); \
    } \
    virtual std::string GetProtocolDisplayName() override \
    { return GetProtocolName(); }
 
#define AddDecoderClass(T) Filter::DoAddDecoderClass(T::GetProtocolName(), T::CreateInstance)
 
#endif