Lmn resampling

aux/inc/WireCellAux/FrameTools.h

@@ -65,7 +65,7 @@ namespace WireCell {
                   channel_list::iterator ch_begin,
                   channel_list::iterator ch_end, int tbin = 0);
 
-        /// Full feature fill.  Fill array with change from traces as
+        /// Full feature fill.  Fill array with charge from traces as
         /// above Return ordered, unique list of channel, each channel
         /// corresponding to one row.  Columns will span the
         /// tbin_range().

aux/inc/WireCellAux/Resampler.h

@@ -0,0 +1,78 @@
+/** Resample frames to a new sampling period/frequency.
+
+    This applies the LMN method to resample a frame under the interpretation
+    that the original sampling is that of an instantaneous quantity.  That is,
+    the resampling applies interpolation normalization.  See the LMN paper for
+    details.
+
+    The resampling may either target a resampled period or a resampled size
+    (number of samples/ticks).  When the target is a period, it must satisfy the
+    LMN rationality condition.  Whether the target is a sampling period or the
+    number samples, the other quantity is determined by the method.
+
+    CAVEAT: this will resample all traces assuming they are dense and with not
+    tbin offset.  Though, traces need not be all the same size.
+
+    Any frame tags, trace tags or trace summaries are carried forward to the
+    output as-given. 
+ */
+
+#ifndef WIRECELLAUX_RESAMPLER
+#define WIRECELLAUX_RESAMPLER
+
+#include "WireCellIface/IFrameFilter.h"
+#include "WireCellIface/IConfigurable.h"
+#include "WireCellAux/Logger.h"
+#include "WireCellIface/IDFT.h"
+
+namespace WireCell::Aux {
+
+    class Resampler : public Aux::Logger,
+                      public IFrameFilter, public IConfigurable {
+      public:
+        Resampler();
+        virtual ~Resampler();
+
+        virtual void configure(const WireCell::Configuration& config);
+        virtual WireCell::Configuration default_configuration() const;
+
+        virtual bool operator()(const input_pointer& inframe, output_pointer& outframe);
+
+      private:
+
+        // Configure: period
+        //
+        // The target sampling period.
+        double m_period{0};
+
+        // Configure: dft
+        //
+        // Name of the DFT component
+        IDFT::pointer m_dft;
+
+        // Configure: time_pad
+        //
+        // Name for a strategy for padding in the time domain when needed to
+        // reach rationality condition.
+        //
+        // - zero :: pad with zero values (default)
+        // - first :: pad with first time sample value
+        // - last :: pad with last time sample value
+        // - median :: pad with median value
+        // - linear :: pad with linear between first and last sample values
+        std::string m_time_pad{"zero"};
+
+        // Configure: time_sizing
+        //
+        // Name a strategy for truncating the final time domain waveform
+        //
+        // - duration :: approximately retain duration (deafult).
+        // - padded :: include rationality condition padding, duration may change.
+        std::string m_time_sizing{"duration"};
+
+        size_t m_count{0};
+
+    };
+
+}
+#endif

aux/src/FrameTools.cxx

@@ -229,8 +229,8 @@ Aux::channel_list Aux::fill(Array::array_xxf& array,
     auto tbinmm = Aux::tbin_range(traces);
     const size_t ncols = tbinmm.second - tbinmm.first;
     const size_t nrows = std::distance(ch.begin(), ch.end());
-    Array::array_xxf arr = Array::array_xxf::Zero(nrows, ncols);
-    Aux::fill(arr, traces, ch.begin(), ch.end(), tbinmm.first);
+    array = Array::array_xxf::Zero(nrows, ncols);
+    Aux::fill(array, traces, ch.begin(), ch.end(), tbinmm.first);
     return ch;
 }
 

aux/src/Resampler.cxx

@@ -0,0 +1,156 @@
+#include "WireCellAux/Resampler.h"
+#include "WireCellUtil/LMN.h"
+#include "WireCellUtil/Units.h"
+
+#include "WireCellAux/SimpleFrame.h"
+#include "WireCellAux/SimpleTrace.h"
+#include "WireCellAux/FrameTools.h"
+#include "WireCellAux/DftTools.h"
+
+#include "WireCellUtil/NamedFactory.h"
+
+#include <map>
+#include <unordered_set>
+
+WIRECELL_FACTORY(Resampler, WireCell::Aux::Resampler,
+                 WireCell::INamed,
+                 WireCell::IFrameFilter, WireCell::IConfigurable)
+
+using namespace std;
+using namespace WireCell;
+using WireCell::Aux::DftTools::fwd_r2c;
+using WireCell::Aux::DftTools::inv_c2r;
+using namespace WireCell::Aux;
+using WireCell::Aux::SimpleTrace;
+using WireCell::Aux::SimpleFrame;
+
+Aux::Resampler::Resampler()
+    : Aux::Logger("Resampler", "aux")
+{
+}
+
+Aux::Resampler::~Resampler() {}
+
+WireCell::Configuration Aux::Resampler::default_configuration() const
+{
+    Configuration cfg;
+
+    cfg["period"] = m_period;
+    cfg["dft"] = "FftwDFT";
+    cfg["time_pad"] = m_time_pad;
+    cfg["time_sizing"] = m_time_sizing;
+    return cfg;
+}
+
+void Aux::Resampler::configure(const WireCell::Configuration& cfg)
+{
+    m_period = get(cfg, "period", m_period);
+    m_time_pad = get(cfg, "time_pad", m_time_pad);
+    m_time_sizing = get(cfg, "time_sizing", m_time_sizing);
+    auto dft_tn = get<std::string>(cfg, "dft", "FftwDFT");
+    m_dft = Factory::find_tn<IDFT>(dft_tn);
+    log->debug("resample period={}, time_pad={} time_sizing={}",
+               m_period, m_time_pad, m_time_sizing);
+}
+
+
+
+bool Aux::Resampler::operator()(const input_pointer& inframe, output_pointer& outframe)
+{
+    outframe = nullptr;
+    if (!inframe) {
+        log->debug("EOS at call={}", m_count);
+        ++m_count;
+        return true;
+    }
+
+    const double Ts = inframe->tick();
+    const double Tr = m_period;
+
+    if (Ts == Tr) {
+        log->warn("same periods: Ts={} Tr={} at call={} is probably not what you want", Ts, Tr, m_count);
+        outframe = inframe;
+        ++m_count;
+        return true;
+    }
+
+    const size_t Nrat = LMN::rational(Ts, Tr);
+
+    std::unordered_map< std::string, IFrame::trace_list_t> tag_indicies;
+
+    ITrace::vector out_traces;
+    for (const auto& trace : *inframe->traces()) {
+
+        if (trace->tbin()) {
+            raise<ValueError>("currently no support for nonzero tbin (fixme)");
+        }
+
+        auto wave = trace->charge();
+
+        const size_t Ns_orig = wave.size();
+        const size_t Ns_pad = LMN::nbigger(Ns_orig, Nrat);
+        const double duration = Ns_pad * Ts;
+        const size_t Nr = duration / Tr; // check for error?
+
+        wave = LMN::resize(wave, Ns_pad);
+        if (Ns_pad > Ns_orig) {
+            const float first = wave[0];
+            const float last = wave[Ns_orig-1];
+
+            auto start = wave.begin()+Ns_orig;
+
+            if (m_time_pad == "linear") {
+                LMN::fill_linear(start, wave.end(), last, first);
+            }
+            // else if (m_time_pad == "cosine") {
+            //     LMN::fill_cosine(start, wave.end(), last, first);
+            // }
+            else {
+                float pad = 0;      // "zero"
+                if (m_time_pad == "first") {
+                    pad = first;
+                }
+                else if (m_time_pad == "last") {
+                    pad = last;
+                }
+                LMN::fill_constant(start, wave.end(), pad);
+            }
+        }
+
+        auto spec = fwd_r2c(m_dft, wave);
+        spec = LMN::resample(spec, Nr);
+        wave = inv_c2r(m_dft, spec);
+
+        // Interpolation interpretation.
+        double norm = wave.size() / (double)Ns_pad; 
+        Waveform::scale(wave, norm);
+
+        if (m_time_sizing == "duration") {
+            const size_t Nr_unpadded = (int)(Ns_orig * (Ts / Tr));
+            wave = LMN::resize(wave, Nr_unpadded);
+        }
+        // else: "padded" is as-is
+
+        if (out_traces.empty()) {
+            log->debug("first ch={} Ts={} Ns={} Ns_pad={} Nrat={} Tr={} Nr={} Nout={} padding:{}",
+                       trace->channel(), Ts, Ns_orig, Ns_pad, Nrat, Tr, Nr, wave.size(), m_time_pad);
+        }
+        out_traces.push_back(std::make_shared<SimpleTrace>(trace->channel(), 0, wave));
+    }
+
+    auto sf = std::make_shared<SimpleFrame>(inframe->ident(), inframe->time(), out_traces, Tr);
+
+    for (const auto& frame_tag : inframe->frame_tags()) {
+        sf->tag_frame(frame_tag);
+    }
+    for (const auto& trace_tag : inframe->trace_tags()) {
+        auto tl = inframe->tagged_traces(trace_tag);
+        auto ts = inframe->trace_summary(trace_tag);
+        sf->tag_traces(trace_tag, tl, ts);
+    }
+    outframe = sf;
+    log->debug("resample {} traces at call={}", out_traces.size(), m_count);
+
+    ++m_count;
+    return true;
+}

cfg/layers/mids/pdsp/api.jsonnet

@@ -8,21 +8,23 @@ local img = import "api/img.jsonnet";
 // Create a mid API object.  No options supported.
 function(services, params, options={}) {
 
+    local pars = std.mergePatch(params, std.get(options, "params", {})),
+
     anodes :: function()
-        low.anodes(params.geometry.drifts, params.geometry.wires_file),
+        low.anodes(pars.geometry.drifts, pars.geometry.wires_file),
 
     drifter :: function(name="")
         low.drifter(services.random,
-                    low.util.driftsToXregions(params.geometry.drifts),
-                    params.lar, name=name),
+                    low.util.driftsToXregions(pars.geometry.drifts),
+                    pars.lar, name=name),
 
     // track_depos, signal, noise, digitizer
-    sim :: sim(services, params),
+    sim :: sim(services, pars, options),
 
     // nf, sp, dnnroi
-    sigproc :: sigproc(services, params, options),
+    sigproc :: sigproc(services, pars, options),
 
-    img :: img(services, params),
+    img :: img(services, pars),
 }
 
 

cfg/layers/mids/pdsp/api/nf.jsonnet

@@ -8,22 +8,22 @@ local chndbs = {
     "perfect": import "chndb-perfect.jsonnet",
 };
 local filters = import "noise-filters.jsonnet";
-local frs = import "frs.jsonnet";
+local resps = import "resps.jsonnet";
 
-function(services, params, options) function(anode)
-    local pars = std.mergePatch(params, std.get(options, "params", {}));
+
+function(services, params, options={}) function(anode)
 
     local ident = low.util.idents(anode);
 
     // nf uses same fr as sp
-    local fr = frs(pars).nf;
+    local fr = resps(params, options).nf.fr;
 
-    local chndb = chndbs[pars.nf.chndb](anode, pars.nf.binning,
+    local chndb = chndbs[params.nf.chndb](anode, params.nf.binning,
                                           fr, services.dft);
 
     // filter look up
-    local flu = filters(services, pars, anode, chndb);
-    local fls = pars.nf.filters;
+    local flu = filters(services, params, anode, chndb);
+    local fls = params.nf.filters;
 
     local fobj = {
         channel: [flu[nam] for nam in fls.channel],

cfg/layers/mids/pdsp/api/resps.jsonnet

@@ -1,4 +1,4 @@
-function(params) {
+function(params, options = {}) {
     sim: {
         er: {
             type: "ColdElecResponse",
@@ -8,7 +8,7 @@ function(params) {
         fr: {
             type: "FieldResponse",
             name: "sim",
-            data: { filename: params.ductor.field_file }
+            data: { filename: std.get(options, "fields", params.ductor.field_file) }
         },
         rc: {
             type: 'RCResponse',
@@ -19,7 +19,7 @@ function(params) {
         fr: {
             type: "FieldResponse",
             name: "nf",
-            data: { filename: params.nf.field_file }
+            data: { filename: std.get(options, "fields", params.nf.field_file) }
         },
     },
     sp: {
@@ -31,7 +31,7 @@ function(params) {
         fr: {
             type: "FieldResponse",
             name: "sp",
-            data: { filename: params.sp.field_file }
+            data: { filename: std.get(options, "fields", params.sp.field_file) }
         },
     },
 

cfg/layers/mids/pdsp/api/sim.jsonnet

@@ -18,7 +18,7 @@ function(services, params, options={}) {
     // Signal binning may be extended from nominal.
     local sig_binning = params.ductor.binning,
 
-    local res = resps(params).sim,
+    local res = resps(params, options).sim,
 
     // some have more than one
     local short_responses = [ res.er, ],

cfg/layers/mids/pdsp/api/sp.jsonnet

@@ -11,15 +11,16 @@ local resps = import "resps.jsonnet";
 
 // Allow an optional argument "sparse" as this is really an end-user
 // decision.  Higher layers may expose this option to the TLA.
-function(services, params, options={}) function(anode)
-    local pars = std.mergePatch(params, std.get(options, "params", {}));
-    local opts = {sparse:true} + options;
+function(services, params, options = {}) function(anode)
+
+    local sparse = std.get(options, 'sparse', true);
     local ident = low.util.idents(anode);
-    local resolution = pars.digi.resolution;
-    local fullscale = pars.digi.fullscale[1] - pars.digi.fullscale[0];
+
+    local resolution = params.digi.resolution;
+    local fullscale = params.digi.fullscale[1] - params.digi.fullscale[0];
     local ADC_mV_ratio = ((1 << resolution) - 1 ) / fullscale;
 
-    local res = resps(params).sp;
+    local res = resps(params, options).sp;
 
     low.pg.pnode({
         type: 'OmnibusSigProc',
@@ -75,6 +76,6 @@ function(services, params, options={}) function(anode)
             mp2_roi_tag: 'mp2_roi' + ident,
 
             isWrapped: false,
-            sparse : opts.sparse,
+            sparse : sparse,
         },
     }, nin=1, nout=1, uses=[anode, services.dft, res.fr, res.er] + spfilt)