From c4eeb4fb43af18eb6bdffcda8d133db5554ab58b Mon Sep 17 00:00:00 2001 From: Brett Viren Date: Fri, 19 Jan 2024 15:57:19 -0500 Subject: [PATCH] Initial support for resampling FRs. Needs more checking --- wirecell/resp/__main__.py | 54 +++---- wirecell/resp/resample.py | 52 ++++++ wirecell/util/__main__.py | 58 +++++++ wirecell/util/functions.py | 3 +- wirecell/util/lmn.py | 323 +++++++++++++++++++++++++++++++++++++ 5 files changed, 456 insertions(+), 34 deletions(-) create mode 100644 wirecell/resp/resample.py create mode 100644 wirecell/util/lmn.py diff --git a/wirecell/resp/__main__.py b/wirecell/resp/__main__.py index e55c9ce..842b5c1 100644 --- a/wirecell/resp/__main__.py +++ b/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()) diff --git a/wirecell/resp/resample.py b/wirecell/resp/resample.py new file mode 100644 index 0000000..e957fc3 --- /dev/null +++ b/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) diff --git a/wirecell/util/__main__.py b/wirecell/util/__main__.py index e932e04..75d1f0c 100644 --- a/wirecell/util/__main__.py +++ b/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") diff --git a/wirecell/util/functions.py b/wirecell/util/functions.py index 8e39297..f052f67 100644 --- a/wirecell/util/functions.py +++ b/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): diff --git a/wirecell/util/lmn.py b/wirecell/util/lmn.py new file mode 100644 index 0000000..e2a3a98 --- /dev/null +++ b/wirecell/util/lmn.py @@ -0,0 +1,323 @@ +#!/usr/bin/env python +''' +The LMN resampling method from the paper . +''' + +import numpy +import math +from numpy import pi +import dataclasses +import matplotlib.pyplot as plt + +@dataclasses.dataclass +class Sampling: + ''' + Characterize a discrete sampling of a signal as a finite sequence. + ''' + + T: float + ''' + The sampling time period + ''' + + N: int + ''' + The number of samples + ''' + + @property + def F(self): + 'The maximum sampling frequency' + return 1.0/self.T + + @property + def duration(self): + 'The time duration of the signal' + return self.N*self.T + + @property + def dF(self): + 'Bandwidth of a frequency sample' + return 1.0/(self.duration) + + @property + def nyquist(self): + 'The Nyquist frequency' + return self.F/2.0 + + @property + def times(self): + 'Sequence of sample times' + return numpy.linspace(0, self.duration, self.N, endpoint=False) + + @property + def freqs(self): + 'Sequence of sample frequencies over the "Nyquist-centered" range.' + return numpy.linspace(0, self.F, self.N, endpoint=False) + + @property + def freqs_zc(self): + 'Sequence of sample frequencies over the "zero-centered" range.' + if self.N%2: # odd + # -F/2-->|[...n_half...][0][...n_half...]|<--F/2 + return numpy.linspace(-self.nyquist, self.nyquist, self.N) + else: # even + return numpy.linspace(-self.nyquist+self.dF, self.nyquist, self.N) + + def resampling(self, N): + ''' + Return a new Sampling with same duration and with N samples. + ''' + if not isinstance(N, int): + raise ValueError(f'N must be integer, got "{type(N)}"') + T = self.duration / N + return Sampling(T, N) + + def __str__(self): + return f'N={self.N} T={self.T}' + +@dataclasses.dataclass +class Signal: + ''' + A sampled signal in time and frequency domain sequence representation. + ''' + + sampling: Sampling + ''' + The time and frequency domain sampling. + ''' + + wave: numpy.ndarray + ''' + Sequence of time domain samples + ''' + + spec: numpy.ndarray + ''' + Sequence of frequency domain samples + ''' + + name: str = "" + ''' + A descriptor of this signal. + ''' + + def __init__(self, sampling, wave=None, spec=None, name=None): + ''' + Construct a signal from a wave, a spectrum or both. + + If one is missing the other is calcualted + ''' + self.sampling = sampling + if wave is None and spec is None: + raise ValueError("must provide at least one of wave or spec") + self.wave = wave + self.spec = spec + self.name = name or str(sampling) + + if wave is None: + self.wave = numpy.real(numpy.fft.ifft(self.spec)) + if spec is None: + self.spec = numpy.fft.fft(self.wave) + + def __str__(self): + return f'{self.name} {sampling}' + + @property + def time_energy(self): + ''' + The amount of Parseval energy in the time domain. + ''' + return numpy.sum( numpy.abs( self.wave * numpy.conj(self.wave) ) ) + + @property + def freq_energy(self): + ''' + The amount of Parseval energy in the frequency domain. + ''' + return numpy.sum( numpy.abs( self.spec * numpy.conj(self.spec) ) ) / self.sampling.N + + +def egcd(a, b, eps=1e-6): + '''Greated common divisor of floating point values a and b. + + Uses the extended Euclidean algorithm. + + A non-zero error is required if values that are not integer valued. + + ''' + def step(a, b): + if a < eps: + return (b, 0, 1) + else: + g, x, y = step(b % a, a) + return (g, y - (b // a) * x, x) + return step(a,b)[0] + + + +def resize_duration(sam, duration, eps=1e-6): + ''' + Return a new sampling with the given duration. + ''' + Ns = sam.N * duration/sam.duration + if abs(Ns - round(Ns)) > eps: + raise ValueError(f'Resizing from duration {sam.duration} to {duration} requires non-integer number of samples of {sam.N} to {Ns}') + Ns = round(newNs) + return Sampling(sam.T, Ns) + + +def resize(sig, duration, pad=0, name=None): + '''Return a new signal of the given duration. + + Note, if duration is not an integer multiple of sig's T, the number + resulting number of samples will be rounded up. Check the .duration of the + result if this matter. + + The signal wave will be truncation or extended with values of pad to fit the + new duration. + ''' + ss = sig.sampling + Nr = math.ceil(ss.N * duration/ss.duration) + + wave = sig.wave.copy() + wave.resize(Nr) # in place + if Nr > ss.N: + wave[ss.N:] = pad + + rs = Sampling(ss.T, Nr) + return Signal(rs, wave=wave, name = name or sig.name) + + +def rational_size(Ts, Tr, eps=1e-6): + ''' + Return a minimum size allowing LMN resampling from period Ts to Tr to be rational. + ''' + dT = Ts - Tr + + n = dT/egcd(Tr, dT, eps=eps) + rn = round(n) + + err = abs(n - rn) + if err > eps: + raise ValueError(f'no GCD for {Tr=}, {Ts=} within error: {err} > {eps}') + + Ns = rn * Tr / dT + rNs = round(Ns) + err = abs(rNs - Ns) + if err > eps: + raise ValueError(f'rationality not met for {Tr=}, {Ts=} within error: {err} > {eps}') + + return rNs + + +def rational(sig, Tr, eps=1e-6): + ''' + Return a new signal like sig but maybe extended to have rational size. + ''' + Ts = sig.sampling.T + nrat = rational_size(Ts, Tr, eps) + Ns = sig.sampling.N + nrag = Ns % nrat + if not nrag: + return sig + + cur = sig.wave.copy() + npad = nrat - nrag + Ns += npad + cur.resize(Ns) + + ss = Sampling(Ts, cur.size) + sig = Signal(ss, wave=cur) + return sig + + +def condition(signal, Tr, eps=1e-6, name=None): + '''Return a signal that is conditioned to be resampled to period Tr. + + The wave will be extended to accomodate a half-cosine with a frequency that + of the peak of the input spectrum. It is further padded so that the entire + waveform is a multiple of the minimimum size for optimal LMN resampling of + the signal to period Tr. Any additional samples are given the amplitude + matching the first sample. + ''' + Ts = signal.sampling.T + Ns = signal.sampling.N + + first = signal.wave[0] + last = signal.wave[-1] + ipeak = numpy.argmax(numpy.abs(signal.spec)) + fpeak = signal.sampling.freqs[ipeak] + + # Number of time bins needed for a half-cosine at fpeak. + nsmooth = int(numpy.ceil(1/(2*Ts*fpeak))) + #print(f'{fpeak=} {Ts=} {nsmooth=}') + + # need at least this much to fit in smoothing half-cosine. + nmin_needed = Ns + nsmooth + + # Final size must be a multiple of this. + nmin = rational_size(Ts, Tr) + + # total padded size to reach next multiple of nmin + npadded = int(numpy.ceil(nmin_needed/nmin)*nmin) + # print(f'{nmin=} {nmin_needed=} {npadded=}') + + wave = numpy.zeros(npadded, dtype=signal.wave.dtype) + wave[:Ns] = signal.wave + wave[Ns:] = first + + # half-cosine parameters + amp = 0.5*(last - first) + bl = 0.5*(first + last) + # the times for making smoothing half-cosine + smooth_time = numpy.linspace(0, nsmooth*Ts, nsmooth, endpoint=False) + smoother = bl + amp*numpy.cos(2*pi*fpeak*smooth_time) + wave[Ns: Ns+nsmooth] = smoother + return Signal(Sampling(T=Ts, N=npadded), wave=wave, name=name) + +def resample(signal, Nr, name=None): + ''' + Return a new signal of same duration that is resampled to have number of samples Nr. + ''' + Ns = signal.sampling.N + resampling = signal.sampling.resampling(Nr) + Tr = resampling.T + + S = signal.spec + R = numpy.zeros(Nr, dtype=S.dtype) + + # Number of unique bins in teh half spectrum ignoring the "DC bin" and + # ignoring the possible Nyquist bin. + if Ns%2: # odd, no nyquist bin + Ns_half = (Ns-1)//2 + else: # even, exclude nyquist bin + Ns_half = (Ns-2)//2 + if Nr%2: # odd, no nyquist bin + Nr_half = (Nr-1)//2 + else: # even, exclude nyquist bin + Nr_half = (Nr-2)//2 + + # Get "half spectrum" size for the non-zero part of resampled spectrum. + if Nr > Ns: # upsample + n_half = Ns_half + else: # downsample + n_half = Nr_half + + # copy DC term and "positive" aka "low-frequency" half spectrum. + R[:1+n_half] = S[:1+n_half] + # copy the "negative" aka "high-frequency" half spectrum + R[-n_half:] = S[-n_half:] + + #print(f'{R.shape=} {Nr_half=} {S.shape=} {Ns_half=} {n_half=}') + + # Deal with Nyquist bins + if Nr > Ns and not Ns%2: # upsampling from spectrum with Nyquist bin + val = S[Ns//2] + R[1+n_half+1] = 0.5*val + R[-n_half-1] = 0.5*val + if Nr < Ns and not Nr%2: # downsampling to spectrum with Nyquist bin + R[1+n_half+1] = abs(S[1+n_half+1]) + + return Signal(resampling, spec = R, name=name) +