Initial support for resampling FRs.

Needs more checking

wirecell/resp/__main__.py

@@ -6,6 +6,7 @@
 import click
 from wirecell.util.fileio import load as source_loader
 from wirecell import units
+from wirecell.util.functions import unitify
 import numpy
 
 cmddef = dict(context_settings = dict(help_option_names=['-h', '--help']))
@@ -58,55 +59,42 @@ def gf_info(dataset):
     dsdict_dump(ds)
 
 
-@cli.command("field-response-transform")
-@click.option("-g", "--gain", default=None, type=str,
-                  help="Set gain with units, eg '14*mV/fC'")
-@click.option("-s", "--shaping", default=None, type=str,
-                  help="Set shaping time with units, eg '2*us'")
-@click.option("-n", "--number", default=None, type=int,
-              help="Resample the field response to have this number of ticks")
+@cli.command("resample")
 @click.option("-t", "--tick", default=None, type=str,
               help="Resample the field response to have this sample period with units, eg '64*ns'")
+@click.option("-P", "--period", default=None, type=str,
+              help="Override the sampling period given in the frfile, eg '100*ns'")
 @click.option("-o", "--output", default="/dev/stdout",
               help="File in which to write the result")
-@click.option("-s", "--schema", default="json.gz", 
-              help="The format to assume given as an extension")
-@click.option("-z", "--compression", default=True,
-              help="Apply compression to the output")
+@click.option("-e", "--error", default=1e-6,
+              help="Precision by which integer and rationality conditions are judged")
 @click.argument("frfile")
-def field_response_transform(gain, shaping, number, tick, output, schema, frfile):
-    '''Apply a transformation to a field response (FR) file.
+def resample(tick, period, output, error, frfile):
+    '''Resample the FR.
 
-    This may be used to transform FR file format or apply a resampling in time.
+    The initial sampling period Ts (fr.period or --period if given) and the
+    resampled period Tr (--tick) must satisfy the LMN rationality condition.
 
-    If both gain and shaping are given then convolve with Cold Electronics
-    response.
+    The total duration of the resampled responses may change.
 
-    If number or tick but not both are given then the duration of the input FR
-    signals is held fixed in the resampling.  If both are given, the duration
-    will made "number*tick".  Resampling must be "lmn-rational".
+    See also:
 
-    See also: wirecell-sigproc fr2npz
+    - wirecell-util lmn 
+    - wirecell-sigproc fr2npz 
 
     '''
+
     import wirecell.sigproc.response.persist as per
-    import wirecell.sigproc.response.arrays as arrs
+    import wirecell.resp.resample as res
 
-    fr = per.load(frfile)
-    if gain and shaping:
-        gain = unitify(gain)
-        shaping = unitify(shaping)
-        dat = arrs.fr2arrays(fr, gain, shaping)
-    else:
-        dat = arrs.fr2arrays(fr)
+    tick = unitify(tick)
+    period = unitify(period)
 
-    if number or tick:
-        tick = unitify(tick) if tick else None
-        dat = arrs.resample(dat, tick, number)
+    fr = per.load(frfile)
+    fr = res.resample(fr, tick, period)
 
-    per.dump(output, dat, 
+    per.dump(output, fr)
 
-    numpy.savez(npz_file, **dat)
 
 def main():
     cli(obj=dict())

wirecell/resp/resample.py

@@ -0,0 +1,52 @@
+#!/usr/bin/env python
+'''
+Functions to resample FR
+'''
+import numpy
+from wirecell import units
+from wirecell.util import lmn
+from wirecell.sigproc.response.schema import FieldResponse, PlaneResponse, PathResponse
+
+
+def resample_one(pr, Ts, Tr, eps=1e-6):
+    '''
+    Resample the path response pr period Ts to Tr.
+    '''
+    sig = lmn.Signal(lmn.Sampling(Ts, pr.current.size), wave = numpy.array(pr.current))
+
+    rat = lmn.rational(sig, Tr, eps)
+
+    Nr = rat.sampling.duration / Tr
+    if abs(Nr - round(Nr)) > eps:
+        raise LogicError('rational resize is not rational')
+    Nr = round(Nr)
+
+    res = lmn.resample(rat, Nr)
+
+    rez = lmn.resize(res, sig.sampling.duration)
+
+    pr = PathResponse(rez.wave, pr.pitchpos, pr.wirepos)
+    return pr
+
+
+def resample(fr, Tr=None, Ts=None, eps=1e-6):
+    '''Return a resampled version of the fr in schema form.
+
+    Either of both of resampled period Tr or size Nr must not be None.  If one
+    is None it will be set so that the total signal duration is retained.  If
+    both are given the resampled signal will have different duration.
+
+    If Ts is given, it overrides the period given in the FR.
+
+    '''
+
+    Ts = Ts or fr.period        # beware, fr.period is sometimes imprecise 
+
+    planes = []
+    for plane in fr.planes:
+        paths = []
+        for pr in plane.paths:
+            pr = resample_one(pr, Ts, Tr, eps=eps)
+            paths.append(pr)
+        planes.append(PlaneResponse(paths, plane.planeid, plane.location, plane.pitch))
+    return FieldResponse(planes, fr.axis, fr.origin, fr.tstart, Tr, fr.speed)

wirecell/util/__main__.py

@@ -18,6 +18,64 @@ def cli(ctx):
     '''
     pass
 
+@cli.command("lmn")
+@click.option("-n", "--sampling-number", default=None, type=int,
+              help="Original number of samples.")
+@click.option("-t", "--sampling-period", default=None, type=str,
+              help="Original sample period,eg '100*ns'")
+@click.option("-T", "--resampling-period", default=None, type=str,
+              help="Resampled sample period, eg '64*ns'")
+@click.option("-e", "--error", default=1e-6,
+              help="Precision by which integer and rationality conditions are judged")
+def cmd_lmn(sampling_number, sampling_period, resampling_period, error):
+    '''Print various LMN parameters for a given resampling.
+
+    '''
+    Ns, Ts, Tr = sampling_number, sampling_period, resampling_period
+
+    if not Ts or not Ns or not Tr:
+        raise click.BadParameter('Must provide all of -n, -t and -T')
+
+    Ts = unitify(Ts)
+    Tr = unitify(Tr)
+
+    from wirecell.util import lmn
+    print(f'initial sampling: {Ns=} {Ts=}')
+
+    nrat = lmn.rational_size(Ts, Tr, error)
+    print(f'rationality factor: {nrat}')
+
+    nrag = Ns % nrat
+    if nrag:
+        npad = nrat - nrag
+        Ns_rational = Ns + npad
+        print(f'rationality padding: {Ns=} += {npad=} -> {Ns_rational=}')
+    else:
+        print(f'rationality met: {Ns=}')
+        Ns_rational = Ns
+
+    Nr_target = Ns_rational*Ts/Tr
+    assert abs(Nr_target - round(Nr_target)) < error # assured by rational_size()
+    Nr_target = round(Nr_target)
+    print(f'resampling target: {Nr_target=} {Tr=}')
+
+    ndiff = Nr_target - Ns_rational
+    print(f'final resampling: {Ns_rational=}, {Ts=} -> {Nr_target=}, {Tr=} diff of {ndiff}')
+
+
+    Nr_wanted = Ns*Ts/Tr
+    Nr = math.ceil(Nr_wanted)
+    if abs(Nr_wanted - Nr) > error:
+        print(f'--> warning: noninteger {Nr_wanted=} for {Tr=}.  Duration will change from {Ns*Ts} to {Nr*Tr} due to rounding.')
+    print(f'requested resampling target: {Nr=} {Tr=}')
+
+    ntrunc = Nr - Nr_target
+    if ntrunc < 0:
+        print(f'truncate resampled: {Nr_target} -> {Nr}, remove {-ntrunc}')
+    if ntrunc > 0:
+        print(f'extend resampled: {Nr_target} -> {Nr}, insert {ntrunc}')
+
+
 
 @cli.command("convert-oneside-wires")
 @click.argument("input-file")

wirecell/util/functions.py

@@ -13,7 +13,8 @@ def unitify(val, unit=""):
     When val is a list, tuple, or dict of unit values, the same
     structure is returned with strings evaluted.
     '''
-
+    if val is None:
+        return None
     if isinstance(val, list):
         return [unitify(one, unit) for one in val]
     if isinstance(val, tuple):