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

Timescale #900

Merged
merged 6 commits into from
May 3, 2023
Merged

Timescale #900

merged 6 commits into from
May 3, 2023

Conversation

FranzForstmayr
Copy link
Collaborator

This PR fixes several issues regarding time transformation.

to be continued

@FranzForstmayr
Copy link
Collaborator Author

Some ongoing work to match the t vector from impulse_response (which looks good) with Frequency.t

@mhuser
Copy link
Collaborator

mhuser commented Apr 25, 2023

Some ongoing work to match the t vector from impulse_response (which looks good) with Frequency.t

Beware pitfall with fft vs rfft. They do not return the same number of samples. I often stumble inside 😉
impulse_response rely on irfft, while time function relies on ifft. This is correct because they do not assume a DC point, but this change the number of samples between the two concepts.

ir = mf.irfft(w.s, n=n)

@FranzForstmayr
Copy link
Collaborator Author

I am aware of the computational difference, but both methods expect the first sample to be the DC component.
In a band pass measurement I'm not sure if the interpretation of time signals make any physical sense without DC extrapolation.
Correct me if I'm wrong.

@FranzForstmayr
Copy link
Collaborator Author

Copied from #860

# Current behaviour
f = rf.Frequency(10, 50, 5)
f.t 
array([-5.0e-11, -2.5e-11,  0.0e+00,  2.5e-11,  5.0e-11]) # 25 ps would imply 40 GSps = 20GHz

# PR 900 now
f = rf.Frequency(10, 50, 5)
f.t 
array([-4.0e-11, -2.0e-11,  0.0e+00,  2.0e-11,  4.0e-11]) # 20 ps would imply 50 GSps = 25 GHz

# I would like to have when PR900 is finished
f = rf.Frequency(10, 50, 5)
f.t 
array([-5.0e-11, -4.0e-11, -3.0e-11, -2.0e-11, -1.0e-11 0.0e+00,  1.0e-11, 2.0e-11, 3.0e-11, 4.0e-11]) # 10 ps would imply 100 GSps = 100 GHz

@mhuser
Copy link
Collaborator

mhuser commented Apr 25, 2023

In a band pass measurement I'm not sure if the interpretation of time signals make any physical sense without DC extrapolation.

According to the Agilent note presented by @Vinc0110 in #757, the bandpass mode without dc still show location and amplitude of discontinuities in impulse response, but not their capacitive or inductive nature. The frequency requirement are less stringent, e.g. no harmonic sweep required.

image

@FranzForstmayr
Copy link
Collaborator Author

Can you calculate this with an ordinary ifft? I think one would need something like STFT to calculatate arbitrary ranges. We do not pass the frequency vector to the ifft currently, so we assume frequency is starting at 0.

@mhuser
Copy link
Collaborator

mhuser commented Apr 25, 2023

Can you calculate this with an ordinary ifft

I think yes. The frequency point still needs to be equally spaced, and thanks to this there is no need to pass the time or frequency axis to the fft, ifft. It can be seen as a time or frequency that is normalized.

In bandpass mode, the ifft time domain results are complex numbers because there is no complex conjugation around a dc center in the input vector; usually, only the magnitude is plotted and it looks terrible, but the discontinuities are supposed to show proper spacing. The frequency sweep is more practical to set around the band of interests if compared to the harmonic sweep that needs to be equally spaced from DC point for lowpass.

Here is another interesting source from Anritsu: Time Domain Measurements Using Vector Network Analyzers pp. 8-10

With a comparison of the methods on a Beatty standard:
image

Prior to it, there is a warning: As noted, Low-pass processing should be used whenever possible as it offers the highest performance and most versatile set of displays

I think I am out of knowledge here because I stick to this advice and always use low-pass mode with dc extrapolation to hopefully understand what I see. I do not truly understand how the bandpass mode works.

In this source, the frequency step for bandpass is (p. 8) df = (f[-1] - f[0]) / N with N + 1 data points.

@mhuser
Copy link
Collaborator

mhuser commented Apr 26, 2023

https://en.wikipedia.org/wiki/Negative_frequency
I still do not get it 100%, but it seems that things are like that:

  • For a complex conjugate Dirac at a given negative and positive same frequency, imaginary sinusoid cancels leaving only real-valued sinusoid, that has the time-domain physical meaning of measured signals
  • Now if there is only a Dirac, e.g. for a positive frequency and zero at negative frequency, there will still be sinusoid but in both real and imaginary domain, with a phase shift between them.

scikit-rf/skrf/network.py

Lines 299 to 301 in c3b5885

'time': mf.ifft,
'time_db': lambda x: mf.complex_2_db(mf.ifft(x)),
'time_mag': lambda x: mf.complex_2_magnitude(mf.ifft(x)),

But in the above code to plot the time domain of s-parameters, which should be equivalent to bandpass, it is like half of the frequency vector will be considered as negative frequencies, which left me confused. It probably works, but how?

@FranzForstmayr
Copy link
Collaborator Author

But in the above code to plot the time domain of s-parameters, which should be equivalent to bandpass, it is like half of the frequency vector will be considered as negative frequencies, which left me confused. It probably works, but how?

I'm also not sure how to interprete the result.

In a band pass measurement I'm not sure if the interpretation of time signals make any physical sense without DC extrapolation.

The current implementation does not care if I pass 0-5GHz, 1-6GHz or 50-55GHz, which seems strange to me.

@mhuser
Copy link
Collaborator

mhuser commented Apr 26, 2023

But in the above code to plot the time domain of s-parameters, which should be equivalent to bandpass, it is like half of the frequency vector will be considered as negative frequencies, which left me confused. It probably works, but how?

I'm also not sure how to interprete the result.

Ok, in this source I think I found half of the answer:
Du to the periodic nature of fft (that's also why we can shift it), it can also be seen as positive-only frequency vector from 0 to Fs. This would be consistant with skrf implementation.
image

But it does not answer our question on how to interpret the bandpass nature of the frequency vector when it start far away 0.

@Vinc0110
Copy link
Collaborator

Hi guys.
I'm loving the discussion, but I currently don't have time to dig into the details. Hopefully I can have a closer look over the weekend.

This whole uncertainty is exactly what irritated me in issue #860. A benchmark result from another tool - known to be correct - would be perfect to compare our implementation against.

@mhuser
Copy link
Collaborator

mhuser commented Apr 27, 2023

This whole uncertainty is exactly what irritated me in issue #860. A benchmark result from another tool - known to be correct - would be perfect to compare our implementation against.

Maybe a brave VNA measurer could measure a length + short and save .s1p as well as bandpass time serie?

In the meantime, I have found that bandpass frequency data seems to be related to downsampling if the fft coefficients are 0 between 0 and a start frequency. The source Oppenheim, Schafer, Discrete-Time Signal Processing is often cited and contains following figure:
image

Unfortunately, I did not found in the book neither how the time vector is reconstructed if the inverse fft is done nor how the time data are supposed to looks like.

@FranzForstmayr
Copy link
Collaborator Author

We could also make a separate Issue/PR for the bandpass stuff and finish this one.
This PR fixes a one off error for the window, the sampling frequency for baseband measured time visualizations and passes the window specification through plot_s_db_time for now.

@mhuser
Copy link
Collaborator

mhuser commented Apr 27, 2023

I think I have got it. I asked a colleague that is much more skilled in signal processing than me and he gave me the interpretation that was later confirmed by the appendix B of Agilent AN 1287-8 Simplified Filter Tuning Using Time Domain

Radio signal transmission - e.g. AM - modulates a carrier frequency with data, which creates sidebands around it. For demodulation, the time-domain data can be multiplied by the carrier frequency, which centers the side bands around 0 Hz and then low-passed to remove the higher-frequencies content.
Bandpass TDR just does that, but because we already are in the frequency domain we just have to pretend that the center frequency of the band is 0 Hz et voilà, demodulation done. Because the spectrum is not symmetric, this will give complex data representing the demodulated I and Q signals. Hence, to locate discontinuities, magnitude should be observed.

In this mode, the center frequency of the frequency
sweep is effectively translated to DC, and the
inverse Fourier transform is applied from minus
one-half of the frequency span to plus one-half of
the span. This is important when looking at a band-
pass filter with a frequency response that is the
same as a low-pass filter response translated up in
frequency to the center of the bandpass filter.
The time-domain transform represents the return
loss as a function of length through the device
under test. For time-domain transforms to be
useful, they must have enough resolution to resolve
the distinguishing characteristics of the network
being measured. In general, the resolution of a
transform is inversely proportional to the frequency
span, although in bandpass mode the resolution
is reduced by half because half the span is for neg-
ative frequencies and half for positive frequencies.

@mhuser
Copy link
Collaborator

mhuser commented Apr 27, 2023

We could also make a separate Issue/PR for the bandpass stuff and finish this one.

Indeed. Maybe should we still wait Monday to let @Vinc0110 have a look if he meet the occasion.

@mhuser
Copy link
Collaborator

mhuser commented Apr 27, 2023

import skrf as rf
import numpy as np

f0 = 10e9
f1 = 50e9
n = 5
f = rf.Frequency(f0, f1, n, unit = 'Hz')

# bandpass frequency and time step
# (according to the app note that the sampling frequency is the band span)
df = (f1 - f0) / (n - 1)
dt = 1 / df

# fft frequency vector
ff1 = np.array([-2, -1, 0, 1, 2]) * df / 5
print(f'{ff1} [Hz]')
ff2 = np.fft.fftshift(np.fft.fftfreq(n, 1/df))
print(f'{ff2} [Hz]')

# fft time vector
tt1 = np.array([-2, -1, 0, 1, 2]) * dt / 5
print(f'{tt1} [s]')
tt2 = np.fft.fftshift(np.fft.fftfreq(n, 1/dt))
print(f'{tt2} [s]')

# output:
[-4.e+09 -2.e+09  0.e+00  2.e+09  4.e+09] [Hz]
[-4.e+09 -2.e+09  0.e+00  2.e+09  4.e+09] [Hz]
[-4.e-11 -2.e-11  0.e+00  2.e-11  4.e-11] [s]
[-4.e-11 -2.e-11  0.e+00  2.e-11  4.e-11] [s]

@FranzForstmayr
Copy link
Collaborator Author

Ok, thank you for this nice paper :)

But then I think that we serve the values for ifft in a wrong order.

From numpy's doc:

The input should be ordered in the same way as is returned by fft, i.e.,
a[0] should contain the zero frequency term,
a[1:n//2] should contain the positive-frequency terms,
a[n//2 + 1:] should contain the negative-frequency terms, in increasing order starting from the most negative frequency.

We calculate the ìfft`like this:

return npy.fft.fftshift(npy.fft.ifft(x, axis=0), axes=0)

where x[0] is our lowest frequency, which is the most negative one. So I think we have to ifftshift the values first.

Copy link
Collaborator

@mhuser mhuser left a comment

Choose a reason for hiding this comment

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

Best if my understanding would be confirmed by someone else.

skrf/frequency.py Outdated Show resolved Hide resolved
skrf/frequency.py Outdated Show resolved Hide resolved
@mhuser
Copy link
Collaborator

mhuser commented Apr 27, 2023

But then I think that we serve the values for ifft in a wrong order.

Oh! I think you are right! Nice spotted!

@mhuser
Copy link
Collaborator

mhuser commented Apr 27, 2023

@FranzForstmayr Can you try this?
(edit: plot without dc extrapolation for bandpass)

import skrf as rf
from skrf.media import Freespace
from matplotlib import pyplot as plt
rf.stylely()

f = rf.Frequency(10, 50, 5, unit = 'GHz')
m = Freespace(f, z0 = 50)
n = m.delay_short(10, unit = 'ps')
n_dc = n.extrapolate_to_dc()
window = 'boxcar'
plt.figure()
plt.title('The 10ps short network')
n_dc.plot_s_time_impulse(window = window, pad = 0, label = 'plot_s_time_impulse')
n_dc.plot_s_time_step(window = window, pad = 0, label = 'plot_s_time_step')
n.plot_s_time_mag(label = 'plot_s_time_mag')

Copy link
Collaborator

@mhuser mhuser left a comment

Choose a reason for hiding this comment

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

@FranzForstmayr
Copy link
Collaborator Author

We already had an issue in the previous version, the first 0 was used as center_to_dc kwarg, that's why two measurements were already plotted. Now the correct measurement is displayed as specified.

@FranzForstmayr
Copy link
Collaborator Author

@FranzForstmayr Can you try this? (edit: plot without dc extrapolation for bandpass)

import skrf as rf
from skrf.media import Freespace
from matplotlib import pyplot as plt
rf.stylely()

f = rf.Frequency(10, 50, 5, unit = 'GHz')
m = Freespace(f, z0 = 50)
n = m.delay_short(10, unit = 'ps')
n_dc = n.extrapolate_to_dc()
window = 'boxcar'
plt.figure()
plt.title('The 10ps short network')
n_dc.plot_s_time_impulse(window = window, pad = 0, label = 'plot_s_time_impulse')
n_dc.plot_s_time_step(window = window, pad = 0, label = 'plot_s_time_step')
n.plot_s_time_mag(label = 'plot_s_time_mag')

This looks good!
test

@mhuser
Copy link
Collaborator

mhuser commented Apr 28, 2023

Great! Nice, the time for the reflexion going back and forth and the resolution divided by two for bandpass are clearly visibles!

@FranzForstmayr FranzForstmayr merged commit ab5ddcb into scikit-rf:master May 3, 2023
10 checks passed
@jhillairet jhillairet mentioned this pull request May 6, 2023
@mhuser mhuser mentioned this pull request May 28, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants