nu-radio
non-causal responses
Describe the issue The RNO-G phased array response as implemented is non-causal
While I haven't checked, I suspect the same is true for many other responses defined by S parameters, unless the S-parameters are derived from time-domain measurements or have been causality-corrected.
Expected behavior
A causal impulse response
Additional context
I think this is caused by 2 issues:
-
S-parameters as measured by network analyzers are usually non-causal (due to numerical errors, etc.). This is also true for the S-parameter files that come from Minicircuits, unfortunately.
-
Even if the s-parameters were causal, linearly interpolating the magnitude/phase would cause it to be non-causal. Causality is guaranteed only when the Kramers-Kronig relations are satisfied, and linear interpolation of the complex response does not preserve that.
I think the solution is either to get a time-domain measurement of the impulse response of the FLOWER, or, at the very least, force causality on resposnse from the S-Parameters. This can be done e.g. by doing complex rational fitting (see e.g. the example here in MATLAB https://www.mathworks.com/help/rf/ref/impulse.html#mw_1b673be6-9aa9-4e4b-bdc7-be829eb756f1 ), which ensures causality, though unfortunately nobody has implemented complex rational fitting (or an equivalent like MATLAB's invfreqs) as part of scipy.signal as far as I can tell. So either someone needs to implement something like it, or we can just use the fitted rational functions from someone doing it manually in MATLAB. Once you have a rational function, you don't need to interpolate either, since your response will be defined at every frequency.
Even if the non-causal part is small, it can have important effects for e.g. template matching or attempts at deconvolution. Plus, I find it extremely embarrassing to plot a non-causal impulse response.
Interesting topic! @fschlueter and I also briefly talked about this in the past and I'm tagging him here in case he is interested in chiming in.
Perhaps this could be a starting point for some experiments?
https://scikit-rf.readthedocs.io/en/latest/api/generated/skrf.vectorFitting.VectorFitting.html
@cozzyd just for my understanding. With "non-causal" in the plot above you refer to the ringing of the trace before t=0? I agree that the linear interpolation I implemented is a very primitive approach. But doing it better, e.g., quadratically, will not make this correct as far as I understand?
The core to the matlab code Cosmin linked has been added to scipy (version >1.15.0) here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.AAA.html#scipy.interpolate.AAA . I'm not seeing any changes
I agree, we should fix this. At very high energies, this also causes artifacts in neutrino simulation as the trigger already fires at a non-causal pre pulse.
@cozzyd did you verify that the effect really comes from your two points above? It should be easy to write a short minimum working example where we just filter a pulse.
A while ago, we also found that cutting the trace to a shorter duration causes low-frequency artifacts which could be mitigated by filtering the beginning and the end of the trace in the time domain with a Hann window (or similar). So maybe the non-causality could also come from another step in the simulation chain(?)
The core to the matlab code Cosmin linked has been added to scipy (version >1.15.0) here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.AAA.html#scipy.interpolate.AAA . I'm not seeing any changes
Hm I am confused now, this still looks non-causal to me with the ringing before the major pulse?
A while ago, we also found that cutting the trace to a shorter duration causes low-frequency artifacts which could be mitigated by filtering the beginning and the end of the trace in the time domain with a Hann window (or similar). So maybe the non-causality could also come from another step in the simulation chain(?)
Wasn't the conclusion back than that this is a "natural effect" which is also observable in data since data is also not windowed?
My two cents (always were) that it comes from the fact that an FFT assumes periodicity and thus is only valid, if you have fully periodic signals. In our case, where we don't, the equivalent of the time-domain treatment (solving coupled differential equations) and an unfolding in Fourier space is not true, hence, you get artifacts, which then need to be suppressed by longer-windows (forcing periodicity to longer time-scales) or windowing. An amp obviously works in the time-domian in real life and not in the frequency domain, so the data should not show that.
Maybe, my point is that we need to identify the real cause and then eliminate it, because (as Cosmin said) it is embaressing to show a non-causal waveform, and (more importantly) it is a real problem for high-energy neutrino simulations. @nilsphysics ran into this problem in his DL reconstruction.
If its artifacts from the FFT, this can be mitigated by time-domain filtering the edges and zero padding. But my feeling is that this is not the main cause of the problem.
The core to the matlab code Cosmin linked has been added to scipy (version >1.15.0) here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.AAA.html#scipy.interpolate.AAA . I'm not seeing any changes
It's not clear to me this does what we want. It seems to be some sort of local approximation for interpolation from what I understand looking quickly at it?
I showed a rational fit of the low-pass at Uppsala last year (using the scikit-rf implementation I linked above). I didn't put the causalized time-domain response on the slides, but it's definitely approximating the S-parameters. I did notice that it's very fiddly and needs some manual tuning to properly converge, so @ryankrebs016 if you want to play with it, be prepared.
Interestingly, the MATLAB RF Toolbox has an SParameter quality tool:
>> sparams = 'ULP-216+_Plus25DegC.s2p'
>> [causality,reciprocity,passivity] = ieee370QualityCheckFrequencyDomain(sparams)
causality =
90.6650
reciprocity =
99.9230
passivity =
100
I'm falling onto the idea that the fitting methods are useful on noisy data where jumps in gain/phase will cause non-causal noise in the impulse response. In this case we may just be doing something weird on this one filter. Is the database detector showing the same issue? I'm only looking at how analog_components.py loads the ULP-216+ response.
After unsuccessfully trying more advanced causal corrections, I checked the spectrums from scikit-rf vs how it's loaded in analog_components.py and saw a difference. Mini-Circuits defines their dB in powers here: https://www.minicircuits.com/appdoc/FILT8-2.html, iirc that means converting dB to a voltage gain we should divide by 20 rather than 10 (scikit-rf and filterresponse.py use 20). Maybe coincidentally then the non-causal ripples nearly disappear. I should be using the same frequency bins as Cosmin's first plot so there shouldn't be different FFT artifacts. Was this the cause? The other Mini-Circuits files with filterresponse.py use 20
wow yeah that looks a lot better... does that mean all S-Parameters are loaded wrong in analog_components?
I am confused, how can a constant factor in the conversion effect the non-physical ripples? However, if you use the rnog_detector.Detector and load the responses from the database, you are using this piece of code defined in the Response class:
if y_ampl_unit == "dB":
gain = 10 ** (y_ampl / 20)
elif y_ampl_unit.lower() == "mag":
gain = y_ampl
else:
raise KeyError
which uses 20.
@ryankrebs016 great! The 10 in the case of the lowpass-filter is definitely incorrect (see bottom of page 6 of the specification for touchstone files).
@fschlueter it's not just an overall constant factor! When you're going from 10^(x/10) to 10^(x/20), you're taking a square root.
Yes I notice that after I wrote it, but we are doing the conversion only in one direction, correct? So we would change the entire normalization?
Yes, it would change the normalization! You can see on Ryan's plots that it does change the normalization a little bit—but, as |S21| ~ 1 in the passband, the change would be small.
okay - sorry for my earlier comments. But this is for the one specific component and the one script we typically do not use? The general problem is not solved. Or is the conclusion that the issue (using somewhat non-causal responses and interpolating them in a way that could introduce further non-causalities) while still present is less severe?
Yes, the underlying problem still exists!
My feeling is that the acausality problem will the strongest in cases where a lot of noise is present in the VNA measurements (as also Ryan notes above). This will be most severe for full-chain measurements where a lot of gain is present in the chain, and the VNA output power has to be correspondingly reduced. In the recent (2024+) measurements, we tried to combat the noise by averaging over many sweeps, but older data is likely to be noisier.
To see what the real impact is, we should probably look directly at the full-chain response from concatenated-S21.
Looks acausal probably due to the "hidden" group delay? FYI I tried fitting with the AAA on this one and it introduced a bunch of acausal noise (so maybe the algorithm doesn't necessarily make clean acausal pulses?).
AAA is a generic rational fitter which can be used for several purposes. It looks like what you linked is probably doing a barycentric rational interpolation, which is not going to enforce causality? We'd want to fit the whole frequency response as one complex laplace-space rational function which will indeed enforce causality.
it looks like perhaps the maximum peak is somehow used to define the 0-point here? This is "fine" as long as it gets accounted for somewhere (cable delays?) but there is also some noise in front of the pulse I think that should at minimum be 0'd out.
Okay, I think I got something working. I'm working off of section 8.2 in "Advanced Signal Integrity for High-Speed Digital Designs." To get nice results I first shifted the pulse to t>0, adding 6ns, to I guess just remove the intrinsic cable delay. Then a new response can be calculated by
new_reponse = np.real(original_reponse) - 1j*np.imag(scipy.signal.Hilbert(np.real(original_response)))
Just to be sure I understand, this is identical to just setting the response to zero for t < 0 (after shifting), correct?
Not claiming it is the thing to do or even smart but I added this function a while ago which determine a "reason on interest" of the response in the time domain and sets everything else to 0.
https://github.com/nu-radio/NuRadioMC/blob/f155b18adcbf54206527497a65d73024bf0df4cc/NuRadioReco/utilities/signal_processing.py#L679
Visibly similar to setting things t<0 equal to zero or smartly windowing but this ensures Kramers-Kronig relations by definition
I'm not totally sure how useful it is since you still need to know what part of the trace is a causal pulse. The intrinsic cable delay in the response is O(0.5ns) whereas I had to shift it by ~6ns