Added a new peak-finding algorithm, fixed saving of filtered frequenc…#997
Added a new peak-finding algorithm, fixed saving of filtered frequenc…#997alisanozdrina merged 3 commits intodevelopfrom
Conversation
fschlueter
left a comment
fschlueter
left a comment
There was a problem hiding this comment.
Seems like a useful addition! I made some comments in-line.
I also have a general question. I am a bit concerned that with a peak_prominance of 4, we are cutting into the thermal noise distribution (the rayleigh distribution has a long tail). Could you make histograms of the amplitudes in a couple of frequency bins, plotted in log and maybe with a rayleigh curve fitted to it (for the raw and the sine subtracted waveforms)? I think with those plots we can clearly see whether that happens. The other question of course if than if that is bad - maybe it is not but at least we should do it internally
| 0.1 to 0.7 GHz is the default for RNO-G, based on bandpass. | ||
| removed_freqs: dict (default: None) | ||
| Dictionary to store filtered noise frequencies for each channel for diagnostics. |
There was a problem hiding this comment.
that is not actually an argument of the function
| Parameters | ||
| ---------- | ||
| ---------- |
There was a problem hiding this comment.
would you mind removing these whitespaces again?
There was a problem hiding this comment.
removed trailing whitespaces
|
|
||
| Returns | ||
| ------- | ||
| -------- |
| algorithm to search for peaks: | ||
| 'simple' search for narrow < 10 MHz peaks which > peak_prominence * RMS withing freq band | ||
| 'sliding' search for narrow < 10 MHz peaks which > peak_prominence * RMS withing 100 MHz sliding window | ||
|
|
|
|
||
| spec_roi = spec[band_mask] | ||
|
|
||
| peak_width_limit = int(10 * 1e-3 / delta_f) #10 MHz |
There was a problem hiding this comment.
you can use 10 * units.MHz instead of 1e-3
There was a problem hiding this comment.
fixed here and one instance later
| spec_roi = spec[band_mask] | ||
| freq_roi = freqs[band_mask] | ||
| # Compute RMS in a sliding window | ||
| window_size = int(50 * 1e-3 / delta_f) #100 MHz |
There was a problem hiding this comment.
comment is wrong, see unit comment
| all_peaks = [] | ||
| for i in range(len(rms_values)): | ||
| start, end = i, i + window_size | ||
| local_spectrum = spec_roi[start:end] |
There was a problem hiding this comment.
can't you use windowed_spectrum here too?
There was a problem hiding this comment.
oh my bad. now fixed
|
I ran sine sub (with new sliding peak searching algo) on ~85K of nuradiomc simulated thermal noise events. Here is the plot. It doesn't cut the tail much, although some bins look a bit different. not sure why though. |
|
Thanks @alisanozdrina those plots are exactly what I head in mind. One question, does the simulation include a realistic bandpass filter, i.e., the channel response? This is with the new peak finder right? Can you also show the old default/way of finding peaks as reference? I would also love to see it on data but I understand if that is to much asked right now. |
fschlueter
added
NuRadioReco
|
For context #885 |
|
Hi @fschlueter, I used Overall, both look good for prominence >= 4; If I set lower prominencies it start cutting the tail and sliding works slightly better. I think it is safe to merge into develop, I'll work on your comments. I didn't have time to run it of forced triggers. If you want test on forced triggers, I would prob do it after new release. 
|
|
@alisanozdrina I've marked this PR as draft for now, once you've addressed Felix' comments you can mark it as "Ready for review" again and/or ping me. |
|
Would be great if we can merge this soon! @alisanozdrina let me elaborate on what I meant with that:
My understanding is that in the simple algorithm you are using a the mean (or median) over a certain frequency range to calculate the threshold for subtracting a sine wave (taking as thresholds a multiple |
…ies for diagnostics.
alisanozdrina
force-pushed
the
channelSinewaveSubtraction_test
branch
from
0777c55 to
2ecfb03
Compare
|
@fschlueter I believe I addressed the comments above. let me know if that's not the case.
Sorry, Do I take it as a green light to merge? Or we still want validation with realistic gain pattern / channel response. |
|
@CamphynR do you think you have the plots in a timely manner? Otherwise I would suggest to merge it asap. |
|
Seeing as the Brussel cluster's storage is completely down at the moment I wouldn't count on it. Merging is ok for me |
|
Give me until this evening hopefully it is resolved by then |
|
I ran with the same settings as Alisa for a limited set of 3 runs of channel 0, station 11, season 2023. These are the results. I might have time later to run it for more forced trigger data but to me these plots do not give a reason to not merge this. The sliding algorithm even brings the distribution closer to a Rayleigh but not sure how relevant that is. 
|
|
@CamphynR, thanks for doing that! yes if it does not take much of your time please run with more data. We try to see whether the tail of the distribution is affected - so more statistics is necessary |
|
Sorry for the delay, these contain a lot more data. You see that only the simple algorithm with peak_promince=3 cuts into the tail (ignore the before in the subtitle that's a typo) |
fschlueter
left a comment
fschlueter
left a comment
There was a problem hiding this comment.
@alisanozdrina and @CamphynR I changed the default method to sliding. Thank you for all your work on this!
4 participants
…ies for diagnostics.
Issue being addressed:
Added an optional peak-selection algorithm to prevent choosing non-CW peaks in the frequency spectrum. New algo is searching for peaks in sliding 100 MHz window instead of the whole spectrum. Tested on Deep CR search dataset burn sample.