Q-factor(s) fitting

[jhillairet]

This is an implementation of the Q-factor (loaded and unloaded) fitting methods proposed by @andrewg22 in #631

This first version is not yet optimized neither fully integrated to scikit-rf (it does not use Network yet)

[jhillairet]

@andrewg22. Currently, the Qfactor() class does not do anything and the code behaves like in the examples, ie with following workflow:

# preliminary: gets f, S and best guesses for Fseed and Qseed

# init Qfactor
Q = rf.Qfactor()

# step1: loaded Q
sv = a_re, a_im, b_re, b_im, QL  = Q.initial_fit(f, S, N, Fseed, Qseed)

# Step 2: Optimised weighted fit
mv, weighting_ratio, number_iterations, RMS_Error = Q.optimise_fit8(f, D, N, Fseed, sv, loop_plan, Tol, quiet)

# unloaded Q-factor and some other useful quantities.
Qo, cal_diam, cal_gamma_V, cal_gamma_T  = Q.Q_unloaded(mv, scaling_factor_A, trmode, quiet)

There is room from improvements, and to make it object-oriented.

Any comments welcome.

[arsenovic]

wow!

[Vinc0110]

Very cool. I've not yet followed the code in much detail, but I will have a closer look later. Is the algorithm doing heavy calculations that require runtime optimizations? If so, I might be able to help with the tricks I learned during the vector fitting improvements.

One very minor thing I noticed: the source code does not contain the word "quality", not a single time. Of course, every rf engineer knows what a Q-factor is, but maybe it's still helpful to mention it once or twice in the docstrings. For example at the beginning:

scikit-rf/skrf/qfactor.py

Line 5 in fd9b29a

Module for fitting Q-factor(s) from S-parameters.

[jhillairet]

Yes, there is room for vectorization improvements, but I would say this is a second step ("premature optimization is..." ;)

What we should discuss first I think is the desired user workflow.

Currently, the fitting is performed by chaining some functions: an initial fit is required before making the improved fitting. Eventually these steps could be combined in a single, but it seems that the precision and the quality of the fit depends much on some user's know-how (like the initial seeds or loop plan)... At this least this is what I understood. @andrewg22 could confirm or not.

[andrewg22]

I've not yet looked at the new code, but it looks as if you are making great progress. I may be able to answer a few points.

The initial and optimised fits can definitely be combined.

I'm sure that some optimisations are possible. Most obviously it could use numpy arrays instead of Python arrays. It's not all that numerically demanding for a typical sweep of 201 points.

It is quite important to have a reasonably accurate initial value for the frequency (Fseed). This can be found from the point of maximum magnitude (transmission) or minimum (reflection or notch resonator). Fseed could be an input parameter with default value None - implying that it will be set from within the function. Qseed is a lot less critical, but it does influence the ability of the optimisation to converge if the VNA sweep is far wider than the resonance. The default value in the test scripts works well.

In some of my work I needed to measure very small shifts in Q-factor, e.g. 0.02 in 400. To allow me to control the iteration accuracy, I provided loop_plan and TOL. The calculation of weighting factors requires a fitted result - so if weightings are required the fitting process must take place at least twice. Loop_plan allows this, but maybe a simpler scheme would be better for SKRF.

[jhillairet]

@andrewg22 and all: I'm refactoring the code and I suggest an API such like that, for comments:

import skrf as rf
ntwk = rf.data.ring_slot

# initalizer. Also perform the initial fit
Qfact = rf.Qfactor(ntwk.s11, res_type='reflection')
print(f' First estimation of Q_L = {Qfact.Q_L}')

# advanced fitting
res = Qfact.fit(method="NLQFIT6", loop_plan='fwfwc')
print(res)

which gives something like:
First estimation of Q_L = [3.39635109]
then

               Q_L: array([3.37972767])
         RMS_Error: array([[0.00095879]])
             f_res: array([8.49194697e+10])
            method: 'NLQFIT6'
 number_iterations: 6
           success: array([[ True]])

what do you think? The optimisation result is returned like a numpy optimization dict.

Should the optimized Q_L be stored in the Qfactor as well? With the same of a different name?

Should the optimized fitting be calculated directly with the initializer (=the inital fit step 1 would then be performed behind the hood)?

[Vinc0110]

Your proposed API seems like the natural approach to me. It actually is very similar to the steps involved in VectorFitting, which might be nice for users doing all sorts of fitting with scikit-rf:

1. Init the fitting algorithm with a Network
2. Call the fitting method
3. Get the results

Only step 3 is different: Qfactor.fit() returns the results directly, whereas VectorFitting.vector_fit() does not return anything but stores the results in its attributes (poles, residues, ...). User can then read them individually or process them with other methods (exporting, plotting, ...). I suggest to do the same for Qfactor.

[jhillairet]

Thanks @Vinc0110.

Another option I also like is for the fit() method to return the fitted version of the object (return self), a solution adopted for example in scikit-learn .

[andrewg22]

@Vinc0110, @jhillairet. These are both good options. I cannot say which would be best.
Should the initialiser have additional parameters? Perhaps the length of uncalibrated line to be de-embedded prior to fitting, and override values for Qseed (default None) and Fseed (default None)?
I don't know if this will influence your thinking, but the fitted object will ideally need some additional methods:
To obtain unloaded Q.
To create the fitted trace - it is very helpful in experiments to upload this into the VNA memory so that it can de compared on its screen. Uncalibrated line (if de-embedded prior to fitting) needs to be re-embedded into the fitted trace.

[jhillairet]

Perhaps the length of uncalibrated line to be de-embedded prior to fitting

For this aspect, scikit-rf has many tools for the simple and advanced deembedding. So my feeling it that the user should provide a Network "finalized", that is, which has been de-embbeded before if necessary using the scikit-rf methods.

and override values for Qseed (default None) and Fseed (default None)?

OK

I don't know if this will influence your thinking, but the fitted object will ideally need some additional methods: To obtain unloaded Q.

Yes definitely. I was wondering if the unloaded Q should be calculated during the call to fit() or with the call to a specific method like Q_unloaded (or similar). What are your thoughts on that?

To create the fitted trace - it is very helpful in experiments to upload this into the VNA memory so that it can de compared on its screen.

OK

Uncalibrated line (if de-embedded prior to fitting) needs to be re-embedded into the fitted trace.

Hum. I need to think about this with respect to what I wrote above...

[andrewg22]

I don't have any strong views on whether to obtain unloaded Q within fit() or by a separate function. It is not always necessary to find unloaded Q, so that would be an argument for making it a separate function .. so that would be a slight argument for making it a separate function.

[jhillairet]

The proposed API is now:

Q = rf.Qfactor(ntwk, res_type='transmission')

# Optimised weighted fit.
# Overwrite Q.Q_L and Q.f_L with fitted results.
res = Q.fit(method="NLQFIT6")

# Unloaded Q-factor and some other useful quantities.
Q0 = Q.Q_unloaded(res, scaling_factor_A)
cal_diam, cal_gamma_V, cal_gamma_T = Q.Q_circle(res, scaling_factor_A)

[jhillairet]

@andrewg22 (but also other of course) , I've added a tutorial:
https://github.com/scikit-rf/scikit-rf/blob/dae54b5e782c5b92f4cb14399c1b14f61960957d/doc/source/tutorials/Q-Factor.ipynb
can you have a look to?

I'm wondering about the API for the fitted trace. Any idea appreciated.

Will also add examples notebook.

[Vinc0110]

That notebook is really nice! I just added a few comments.

[jhillairet]

OK, I think it's ready now.

I've simplified the API. Now the results from the fitting are stored in the class and used per default.

>>> Q = rf.Qfactor(rf.data.ring_slot_meas, res_type='reflection')
>>> opt_res = Q.fit()  # returns the optimized results if needed. Different optimizer can be selected if needed.

>>> print(Q)
Q-factor of Network ring slot measured. (fitted: f_L=85.968GHz, Q_L=3.685)

# resonant frequency and 3dB bandwidth
>>> print(Q.f_L)
[8.59683221e+10]
>>> print(Q.BW)
[2.33318004e+10]

# also exist in scaled to the frequency unit
>>> print(Q.f_L_scaled)  # here in GHz
[85.96832211]
>>> print(Q.BW_scaled)  # here in GHz
[23.33180037]

# passing opt_res is now optional for Unloaded Q, Q-circle, S-parameters and fitted network. 
# If not passed, it uses the one produced during the the last call to fit()
>>> print(Q.Q_unloaded())
[7.79553387]

>>> fitted_ntwk = Q.fitted_network()  

# comparing data and model
>>> rf.data.ring_slot_meas.plot_s_db()
>>> fitted_ntwk.plot_s_db(label='fitted resonant model')
# in polar
>>> rf.data.ring_slot_meas.plot_s_polar()
>>> fitted_ntwk.plot_s_polar(label='fitted resonant model')

[jhillairet]

I've also added an example of application, that is finding relative permittivity and tan(delta) from Q factors fitting, from data provided by Martin Schöön (thanks!).

[jhillairet]

Rendered Doc overview for comments:

https://scikit-rf-jhillairet.readthedocs.io/en/latest/tutorials/Q-Factor.html
https://scikit-rf-jhillairet.readthedocs.io/en/latest/examples/qfactor/Finding_Dk_Df_from_resonance_fitting.html
(weird equation at the end)

[Vinc0110]

Awesome! I also like the new example with multiple resonances. It makes it very clear why is is sometimes necessary to specify the frequency range of interest.

[jhillairet]

I'll merge in the present state. One issue remains on the definition of the bandwidth for reflexion resonator. I've exchanged with a few colleagues on how they define it, but it seems there is not a clear definition. So currently the .BW property returns the value as defined for transmission resonator (resonance freq/Q_L), but I understand this is incorrect.