Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow negative delays in trace_utilities.delay #768

Open
wants to merge 11 commits into
base: develop
Choose a base branch
from
7 changes: 4 additions & 3 deletions NuRadioReco/framework/base_trace.py
Original file line number Diff line number Diff line change
@@ -258,6 +258,7 @@ def apply_time_shift(self, delta_t, silent=False):
"""
if delta_t > .1 * self.get_number_of_samples() / self.get_sampling_rate() and not silent:
logger.warning('Trace is shifted by more than 10% of its length')

spec = self.get_frequency_spectrum()
spec *= np.exp(-2.j * np.pi * delta_t * self.get_frequencies())
self.set_frequency_spectrum(spec, self._sampling_rate)
@@ -299,9 +300,9 @@ def deserialize(self, data_pkl):

def add_to_trace(self, channel):
"""
Adds the trace of another channel to the trace of this channel. The trace is only added within the
Adds the trace of another channel to the trace of this channel. The trace is only added within the
time window of "this" channel.
If this channel is an empty trace with a defined _sampling_rate and _trace_start_time, and a
If this channel is an empty trace with a defined _sampling_rate and _trace_start_time, and a
_time_trace containing zeros, this function can be seen as recording a channel in the specified
readout window.

@@ -354,7 +355,7 @@ def add_to_trace(self, channel):
t_end_readout = tt_readout[i_end_readout]
i_end_channel = n_samples_channel - 1
t_end_channel = t1_channel

# Determine the remaining time between the binning of the two traces and use time shift as interpolation:
residual_time_offset = t_start_channel - t_start_readout
tmp_channel = copy.deepcopy(channel)
12 changes: 6 additions & 6 deletions NuRadioReco/modules/analogToDigitalConverter.py
Original file line number Diff line number Diff line change
@@ -306,12 +306,12 @@ def get_digital_trace(self, station, det, channel,

trace = np.fft.irfft(trace_fft * trigger_filter)

# Random clock offset
delayed_samples = len(trace) - int(np.round(MC_sampling_frequency / adc_sampling_frequency)) - 1
trace = delay_trace(trace, MC_sampling_frequency, adc_time_delay, delayed_samples)

times = times + 1.0 / adc_sampling_frequency
times = times[:len(trace)]
if adc_time_delay:
# Random clock offset
trace, dt_tstart = delay_trace(trace, MC_sampling_frequency, adc_time_delay)
if dt_tstart > 0:
times = times[dt_tstart / MC_sampling_frequency:]
times = times[:len(trace)]

# Upsampling to 5 GHz before downsampling using interpolation.
# We cannot downsample with a Fourier method because we want to keep
88 changes: 60 additions & 28 deletions NuRadioReco/utilities/trace_utilities.py
Original file line number Diff line number Diff line change
@@ -1,10 +1,9 @@
from NuRadioReco.utilities import units, ice, geometryUtilities as geo_utl, fft
import NuRadioReco.framework.base_trace

import numpy as np
import scipy.constants
import scipy.signal
from NuRadioReco.utilities import units
from NuRadioReco.utilities import ice
from NuRadioReco.utilities import geometryUtilities as geo_utl
from NuRadioReco.utilities import fft
import logging
logger = logging.getLogger('NuRadioReco.trace_utilities')

@@ -20,7 +19,7 @@ def get_efield_antenna_factor(station, frequencies, channels, detector, zenith,

Parameters
----------

station: Station
frequencies: array of complex
frequencies of the radio signal for which the antenna response is needed
@@ -195,7 +194,7 @@ def upsampling_fir(trace, original_sampling_frequency, int_factor=2, ntaps=2 **

Parameters
----------

trace: array of floats
Trace to be upsampled
original_sampling_frequency: float
@@ -304,49 +303,82 @@ def apply_butterworth(spectrum, frequencies, passband, order=8):
return filtered_spectrum


def delay_trace(trace, sampling_frequency, time_delay, delayed_samples=None):
def delay_trace(trace, sampling_frequency, time_delay, crop_trace=True):
"""
Delays a trace by transforming it to frequency and multiplying by phases.
Since this method is cyclic, the trace has to be cropped. It only accepts
positive delays, so some samples from the beginning are thrown away and then
some samples from the end so that the total number of samples is equal to

Since this method is cyclic, the delayed trace has to be cropped. A positive
delay means that the trace is shifted to the right, i.e., its delayed. A
negative delay would mean that the trace is shifted to the left. Samples
from the beginning or end are thrown away. Optionally one can crop the
trace (from the end) so that the total number of samples is equal to
the argument delayed samples.

Parameters
----------
trace: array of floats
trace: array of floats or `NuRadioReco.framework.base_trace.BaseTrace`
Array containing the trace
sampling_frequency: float
Sampling rate for the trace
time_delay: float
Time delay used for transforming the trace. Must be positive or 0
delayed_samples: integer or None
Number of samples that the delayed trace must contain
if None: the trace is not cut
crop_trace: bool
If True (default), the trace is cropped to remove samples
what are unphysical after delaying (rolling) the trace.

Returns
-------
delayed_trace: array of floats
The delayed, cropped trace
dt_start: float (optional)
The delta t of the trace start time. Only returned if crop_trace is True.
"""

if time_delay < 0:
msg = 'Time delay must be positive'
raise ValueError(msg)

n_samples = len(trace)

spectrum = fft.time2freq(trace, sampling_frequency)
frequencies = np.fft.rfftfreq(n_samples, 1 / sampling_frequency)
# Do nothing if time_delay is 0
if not time_delay:
if isinstance(trace, NuRadioReco.framework.base_trace.BaseTrace):
if crop_trace:
return trace.get_trace(), 0
else:
return trace.get_trace()
else:
if crop_trace:
return trace, 0
else:
return trace

if isinstance(trace, NuRadioReco.framework.base_trace.BaseTrace):
spectrum = trace.get_frequency_spectrum()
frequencies = trace.get_frequencies()
else:
n_samples = len(trace)
spectrum = fft.time2freq(trace, sampling_frequency)
frequencies = np.fft.rfftfreq(n_samples, 1 / sampling_frequency)

spectrum *= np.exp(-1j * 2 * np.pi * frequencies * time_delay)

delayed_trace = fft.freq2time(spectrum, sampling_frequency)

init_sample = int(time_delay * sampling_frequency) + 1
cycled_samples = int(round(time_delay * sampling_frequency))

if delayed_samples is not None:
delayed_trace = delayed_trace[init_sample:None]
delayed_trace = delayed_trace[:delayed_samples]
if crop_trace:
if time_delay >= 0:
delayed_trace = delayed_trace[cycled_samples:]
dt_start = cycled_samples * sampling_frequency
else:
delayed_trace = delayed_trace[:-cycled_samples]
dt_start = 0

return delayed_trace, dt_start

else:
# Check if unphysical samples contain any signal and if so, throw a warning
if time_delay > 0:
if np.any(np.abs(delayed_trace[:cycled_samples]) > 0.01 * units.microvolt):
logger.warning("The delayed trace has unphysical samples that contain signal. "
"Consider cropping the trace to remove these samples.")
else:
if np.any(np.abs(delayed_trace[-cycled_samples:]) > 0.01 * units.microvolt):
logger.warning("The delayed trace has unphysical samples that contain signal. "
"Consider cropping the trace to remove these samples.")

return delayed_trace
return delayed_trace