Add cw filter notch module #724

CamphynR · 2024-09-19T09:39:10Z

Inclusion of basic cw filter to NuRadioMC which uses notch filters on peaks in the frequency spectrum. The threshold for defining a peak and the width of the notch filter can be modified by parameters in the module begin() function.

This filter would serve as a quick option for people who want to filter CWs in their current work. The notch filter filters out a very narrow band around the frequency peak. To narrow(broaden) the band, increase(decrease) the parameter Q. To be more(less) strict with what to consider a CW, increase(decrease) the threshold parameter.

Note that this filter would not be extremely accurate but can be of use for its speed before the more sophisticated ARA/sine reduction filter is implemented in NuRadio (see /https://github.com/nichol77/libRootFftwWrapper)

fschlueter · 2024-09-19T12:44:09Z

@CamphynR thank you for this module! I have a few general remarks:

I like that you keep the functions which perform the "operations" general such that one can import them from other code and use them there. However, for the sake of the NuRadioReco module you should use trace.get_spectrum (the trace class stores the data internally either in the time or frequency domain depending on the last time it was used and with that you might save one fourier transformation. You can still have a function with does the transformation but this is than not used with in the module (does that make sense?).
please use the sampling rate as stored in the channel object and not hand-set by the module

fschlueter · 2024-09-19T12:44:36Z

NuRadioReco/modules/cwFilterNotch.py

+    return freq[peak_idxs]
+
+
+def filter_cws(trace, Q = 1e3, threshold = 4, fs = 3.2e9 * units.Hz):


please use a different name for Q

fschlueter · 2024-09-19T12:44:49Z

NuRadioReco/modules/cwFilterNotch.py

+    ----------
+    trace : np.ndarray
+        waveform (shape: [24,2048])
+    Q: int, default = 1000


fschlueter · 2024-09-19T12:44:54Z

NuRadioReco/modules/cwFilterNotch.py

+    threshold: int, default = 4
+        threshold for peak definition. A peak is defined as a point in the frequency spectrum
+        that exceeds threshold * rms(real fourier transform)
+    fs: float, default = 3.2e9 Hz


fschlueter · 2024-09-19T12:51:56Z

So for example find_frequency_peaks should take the frequency spectrum but you could add a wrapper to take also the time series.

…between time and frequency domains, also docstring imporvements

CamphynR · 2024-09-23T10:00:29Z

Thanks for the feedback. I tried implementing it in such a way that the code only uses the get_spectrum and get_trace methods to switch between time and frequency domains. Since the notch filter can only be applied directly on the time trace I think it's difficult to leave out any transforms but this should be more streamlined. I also added a trace wrapper for find_frequency_peaks

fschlueter · 2024-09-23T10:06:06Z

NuRadioReco/modules/cwFilterNotch.py

+    ax.legend(loc = legendloc)
+
+
+class cwFilterNotch():


typically we encode in the name on what the module acts on. So for example channelCWNotchFilter

fschlueter · 2024-10-08T07:46:24Z

1

@@ -0,0 +1 @@
+/user/rcamphyn/envs/NuRadioMC: ceph.quota.max_bytes: No such attribute


please remove this file

fschlueter · 2024-10-08T07:48:42Z

NuRadioReco/modules/channelCWNotchFilter.py

+
+"""
+Contains function to filter continuous wave out of the signal
+functions should work on a per event basis to comply with the iteration methods used in readRNOGData


In general all modules work like this. You do not need to mention that here. And I would remove the reference to readRNOGData. This module can also be used for other data/experiments

fschlueter · 2024-10-08T07:49:22Z

NuRadioReco/modules/channelCWNotchFilter.py

+    fs: float, default = 3.2e9 Hz
+        sampling frequency of the RNO-G DAQ, (input should be taking from the channel object)


In line with the comment above I would not refer here to RNO-G. I would probably also remove the default value.

fschlueter · 2024-10-08T07:49:48Z

NuRadioReco/modules/channelCWNotchFilter.py

+    freq_peaks : np.ndarray
+        frequencies at which a peak was found
+    """
+    freq = np.fft.rfftfreq(2048, d = 1/fs)


replace the hardcoded 2048

fschlueter · 2024-10-08T07:53:31Z

NuRadioReco/modules/channelCWNotchFilter.py

+    freq_peaks : np.ndarray
+        frequencies at which a peak was found
+    """
+    freq = np.fft.rfftfreq(2048, d = 1/fs)


Pep8 recommends to not add spaces around the = for keyword arguments: https://peps.python.org/pep-0008/#other-recommendations

This is everywhere in this file. I let you decide if you want to change it

fschlueter · 2024-10-08T07:54:57Z

NuRadioReco/modules/channelCWNotchFilter.py

+    Parameters
+    ----------
+    freq : np.ndarray
+        frequencies of a NuRadio rime trace


fschlueter · 2024-10-08T07:58:51Z

NuRadioReco/modules/channelCWNotchFilter.py

+    """
+    freqs = find_frequency_peaks(freq, spectrum, threshold = threshold)
+
+    if len(freqs) !=0:


space. or simply len(freqs)

fschlueter · 2024-10-08T11:32:15Z

NuRadioReco/modules/channelCWNotchFilter.py

@@ -28,24 +30,24 @@ def find_frequency_peaks_from_trace(trace : np.ndarray, fs : float = 3.2e9 * uni
    freq_peaks : np.ndarray
        frequencies at which a peak was found
    """
-    freq = np.fft.rfftfreq(2048, d = 1/fs)
+    freq = np.fft.rfftfreq(nr_samples, d = 1/fs)


I think len(trace) would be fine here.

fschlueter · 2024-10-08T11:32:59Z

NuRadioReco/modules/channelCWNotchFilter.py

@@ -28,24 +30,24 @@ def find_frequency_peaks_from_trace(trace : np.ndarray, fs : float = 3.2e9 * uni
    freq_peaks : np.ndarray
        frequencies at which a peak was found
    """
-    freq = np.fft.rfftfreq(2048, d = 1/fs)
+    freq = np.fft.rfftfreq(nr_samples, d = 1/fs)
    ft = fft.time2freq(trace, fs)

    freq_peaks = find_frequency_peaks(freq, ft, fs = fs, threshold = threshold)


I think the pep8 rule includes this line too

fschlueter · 2024-10-08T11:51:16Z

@CamphynR thanks for doing this and sorry for my nit-picking! I also made a small commit and refactored a few things, I hope this is okay.

fschlueter · 2024-10-08T11:52:50Z

NuRadioReco/modules/channelCWNotchFilter.py

+def filter_cws(trace : np.ndarray, freq : np.ndarray, spectrum : np.ndarray, fs=3.2e9 * units.Hz, quality_factor=1e3, threshold=4):
+    """
+    Function that applies a notch filter at the frequency peaks of a given time trace
+    using the scipy library
+
+    Parameters
+    ----------
+    trace : np.ndarray
+        Waveform (shape: [24,2048])
+    freq : np.ndarray
+        Frequency of the trace's real fourier transform
+    spectrum:
+        the trace's real fourier transform
+    fs : float, default = 3.2e9 Hz
+        sampling frequency of the RNO-G DAQ
+    quality_factor : int, default = 1000
+        quality factor of the notch filter, defined as the ratio f0/bw, where f0 is the centre frequency
+        and bw the bandwidth of the filter at (f0,-3 dB)
+    threshold : int, default = 4
+        threshold for peak definition. A peak is defined as a point in the frequency spectrum
+        that exceeds threshold * rms(real fourier transform)
+


Is it true that the shape of trace can be [channels, samples]? A few rno-g specific things which could be removed.

You're right the way it's implemented in the class is that it takes a trace from one channel so the shape is actually just [samples]

fschlueter · 2024-10-08T11:53:37Z

Please only answer the question regarding the shape of the trace in the one class. The rest you can ignore and we can approve this

cozzyd

Hi Ruben,

This looks like a useful module!

I'd like to see some clarification on a few things:

Does filtfilt pad the waveform in this case? I don't think that's what we want.
Does the (acausal!) phase response of filtfilt cause any obvious problems with calpulsers? I'm also concerned it might potentially adversely effect e.g. CR template matching.
See if you can get the filtering to work with second-order sections for stability (and this will also let you apply everything once). Since making the comment on the code, I found that https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.sosfiltfilt.html#scipy.signal.sosfiltfilt exists and is probably what you want to use.
I think there can be some simple improvements in the spectral peak finding that could help a lot, like specifying a frequency range for the RMS (otherwise your RMS is basically just the dynamic range of the filter response...). You could also imagine doing things like comparing to a moving-median filtered spectrum or something, though that's probably significantly more expensive. Using some FFT window may also help with finding the right peak.

cozzyd · 2024-10-08T19:39:19Z

NuRadioReco/modules/channelCWNotchFilter.py

+
+    if len(freqs):
+        notch_filters = [signal.iirnotch(freq, quality_factor, fs = fs) for freq in freqs]
+        trace_notched = signal.filtfilt(notch_filters[0][0], notch_filters[0][1], trace)


filtfilt will pad by default, no? (Maybe I'm misremembering)

If it does pad, I see two issues:

you get a longer trace than you put in, which is perhaps surprising

you can get weird transients at the former edge of the waveform, where you'll essentially see the step response of the filter

Hi Cosmin, thank you for the input. I am very much not experienced with signal analysis / filters so I might be asking some obvious questions in what follows.
Filtfit indeed pads the waveform. I looked in the documentation and it applies an "odd padding" (extending the waveform by copying a range of samples at the edge, rotated by 180 degrees). The trace that is returned is the same shape as the input trace so after filter application it cuts the padding. Is it a problem to use padding? In my naive estimation I thought the filter would lead to some artifacts at the edges when not using any padding? Do you mean there would be transients due to the step reponse if the the wavefrom was padded with 0's?

If it does indeed crop the padding applied then that's good and we don't have to worry about that part.

And yes, I'm somewhat worried about the step response of the filter at the edge after padding, since you'll be going from 0 to sometimes a non-zero value that sometimes will be large (i.e. > 3 sigma 0.3% of the time), but I guess as long as we never consider anything on the edge of the waveforms it's probably not a problem?

For now I have set the padding to None (not sure how Scipy then arrives at the same trace length though, but it does), Odd padding is also an option as it basically "extends" the waveform. I do not have a feeling about which option is better.

cozzyd · 2024-10-08T19:49:08Z

NuRadioReco/modules/channelCWNotchFilter.py

+    """
+
+    rms = np.sqrt(np.mean(np.abs(spectrum)**2))
+    peak_idxs = np.where(np.abs(spectrum) > threshold * rms)[0]


Won't the RMS of the trace be dominated by the filter shape? I wonder if it makes sense to limit the frequency range over which the RMS is calculated?

With filter shape do you mean the CW peak?

what I mean is if you take a noiseless neutrino or calpulser with no CW, the RMS will mostly come from the fact that we have a high pass and low pass filters so the spectrum has a shape in it.

cozzyd · 2024-10-09T13:34:26Z

NuRadioReco/modules/channelCWNotchFilter.py

+        notch_filters = [signal.iirnotch(freq, quality_factor, fs = fs) for freq in freqs]
+        trace_notched = signal.filtfilt(notch_filters[0][0], notch_filters[0][1], trace)
+        for notch in notch_filters[1:]:
+            trace_notched = signal.filtfilt(notch[0], notch[1], trace_notched)


I don't know how filtfilt is implemented, but depending on implementation it could make sense to cascade the filters first and only apply once? Though that could be numerically unstable, depending on implementation (I think, what you'd want to do is convert all of the filters into second-order sections and then cascade? I know in MATLAB filtfilt supports filters in second-order section form as input but I don't know how it works in scipy at all, unfortunately. There probably should be away to use second-order sections ('sos') otherwise you run into numerical stability problems at high orders anyway?)

What do you mean with acausal phase response? I was under the impression filtfilt yielded a filter with phase 0 due to it applying once forward and once backwards?

Right, a 0-phase linear filter (other than the trivial one) won't be causal (for example, a brick wall/sinc filter has 0 phase, and has a really acausal response, but in general this should follow from the Kramers-Kronig relations). So this means it could look ugly when plotted and affect time-domain metrics like impulsivity, but on the other hand it has good properties for cross-correlation (it basically just acts as a frequency weighting for cross-correlation, so it mostly just changes the normalization for cross-correlation). Crucially, this won't bias the inter-channel lag measurements (very important if you're applying different filters to different channels), so I guess a 0-phase filter has to be the right choice here!

I wonder if it makes sense to optionally return the applied filter as well so that it also could be applied to e.g. a template used for cross-correlation, or as some sort of correction to an impulsivity measurement? This could be implemented by having an optional kwarg that is either None (default) or a list-like that the filter coefficients could be appended to.

…nce, add logging

fschlueter · 2024-10-16T12:23:19Z

NuRadioReco/modules/channelCWNotchFilter.py

@@ -1,3 +1,6 @@
+import logging
+logger = logging.getLogger(__name__)


can you please check the PR about logging (which gets merged to day) and initialize the logger according to the guidelines introduced there?

cozzyd · 2024-10-16T13:55:08Z

NuRadioReco/modules/channelCWNotchFilter.py

-        trace_notched = signal.filtfilt(notch_filters[0][0], notch_filters[0][1], trace)
-        for notch in notch_filters[1:]:
-            trace_notched = signal.filtfilt(notch[0], notch[1], trace_notched)
+        notch_filters = np.array([signal.iirnotch(freq, quality_factor, fs = fs) for freq in freqs]).reshape(-1, 6)


I don't think this is the right way to convert to SOS form? You probably need to use a filter design method that can return SOS form rather than BA form.

Nevermind, this works because apparently iirnotch always returns a second-order filter!

(I thought somehow that Q might affect the order, but it doesn't...)

…ilter shape

fschlueter · 2024-10-17T08:52:24Z

NuRadioReco/modules/channelCWNotchFilter.py

+    if cache is not None:
+        # Check to avoid cache dictionary overflowing the memory,
+        # set to roughly stay below 6 MB (every filter is 6 floats + freq = 7 floats ~ 56B)
+        if len(cache.keys() < 1e5):


this line seems to be buggy

does the filter depend on the number of samples of the frequency spectrum?

I guess that is handled by the sampling rate which reflects that. I would opt again of removing the default value. Because you might what to upsample the trace before for more accuracy but forget to adjust the sampling rate.

fschlueter · 2024-10-17T09:57:16Z

NuRadioReco/modules/channelCWNotchFilter.py

@@ -91,11 +92,11 @@ def get_filter(freq : int, quality_factor=1e3, fs=3.2e9 * units.Hz, cache=None):
    if cache is not None:
        # Check to avoid cache dictionary overflowing the memory,
        # set to roughly stay below 6 MB (every filter is 6 floats + freq = 7 floats ~ 56B)
-        if len(cache.keys() < 1e5):
+        if len(cache.keys()) < 1e5:


FYI the .keys() is not necessary but it is more verbose.

fschlueter · 2024-11-26T09:29:27Z

@CamphynR can you fix the merge conflict? Also building the documentation fails because of your new module. The failing unit test should be due to another already solved problem and should disappear when you run the tests again.

Cosmetic changes

CamphynR added 2 commits September 19, 2024 11:30

add cw filter which uses notch filters

77dd42c

add description to changelog

275db5f

CamphynR added the enhancement New feature or request label Sep 19, 2024

CamphynR requested review from cozzyd, fschlueter, EnriqueHuesca and sjoerd-bouma September 19, 2024 09:39

fschlueter reviewed Sep 19, 2024

View reviewed changes

attempt to only use trace.get_spectrum and trace.get_trace to switch …

a38751b

…between time and frequency domains, also docstring imporvements

fschlueter reviewed Sep 23, 2024

View reviewed changes

rename to better match NuRadio syntax, bugfix test figure save path

586af75

CamphynR requested a review from fschlueter October 8, 2024 07:31

fschlueter reviewed Oct 8, 2024

View reviewed changes

1 Outdated

@@ -0,0 +1 @@

/user/rcamphyn/envs/NuRadioMC: ceph.quota.max_bytes: No such attribute

Copy link

Collaborator

fschlueter Oct 8, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please remove this file

fschlueter reviewed Oct 8, 2024

View reviewed changes

CamphynR added 2 commits October 8, 2024 13:17

remove file 1

9866711

implement better syntax and fix typos

58ec584

fschlueter reviewed Oct 8, 2024

View reviewed changes

CamphynR and others added 2 commits October 8, 2024 13:39

typo and remove keyword nr_samples

1d40759

Update channelCWNotchFilter.py: Style

211bbbb

fschlueter reviewed Oct 8, 2024

View reviewed changes

change shape in docstring

29d96ca

CamphynR requested a review from fschlueter October 9, 2024 12:15

cozzyd suggested changes Oct 9, 2024

View reviewed changes

Change filter application to apply cascade of second order sections o…

85da819

…nce, add logging

fschlueter reviewed Oct 16, 2024

View reviewed changes

cozzyd reviewed Oct 16, 2024

View reviewed changes

CamphynR added 3 commits October 17, 2024 09:39

make logger compatible with NuRadio logging, add comment explaining f…

51a4334

…ilter shape

add option to save filter coefficients for future reference

31c3806

add caching

3b2d804

fschlueter reviewed Oct 17, 2024

View reviewed changes

CamphynR added 2 commits October 17, 2024 11:45

fix typo

24f6b5c

remove default value for sampling frequency fs

8ee0b00

fschlueter reviewed Oct 17, 2024

View reviewed changes

CamphynR requested a review from cozzyd November 19, 2024 09:21

fschlueter and others added 2 commits November 26, 2024 10:37

Update channelCWNotchFilter.py

ded7d6d

Cosmetic changes

Merge branch 'develop' into add_cwFilterNotch_module

d8a9850

fschlueter approved these changes Nov 26, 2024

View reviewed changes

CamphynR merged commit ff1fd1e into develop Nov 29, 2024
5 checks passed

CamphynR deleted the add_cwFilterNotch_module branch November 29, 2024 12:34

		return freq[peak_idxs]


		def filter_cws(trace, Q = 1e3, threshold = 4, fs = 3.2e9 * units.Hz):

		@@ -0,0 +1 @@
		/user/rcamphyn/envs/NuRadioMC: ceph.quota.max_bytes: No such attribute

		fs: float, default = 3.2e9 Hz
		sampling frequency of the RNO-G DAQ, (input should be taking from the channel object)

		@@ -1,3 +1,6 @@
		import logging
		logger = logging.getLogger(__name__)

Add cw filter notch module #724

Add cw filter notch module #724

Conversation

CamphynR commented Sep 19, 2024

fschlueter commented Sep 19, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fschlueter commented Sep 19, 2024

CamphynR commented Sep 23, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fschlueter commented Oct 8, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fschlueter commented Oct 8, 2024

cozzyd left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

CamphynR Oct 17, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

CamphynR Oct 15, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fschlueter commented Nov 26, 2024

CamphynR Oct 17, 2024 •

edited

Loading

CamphynR Oct 15, 2024 •

edited

Loading