Lorentzian signal fit
What you have for your 2nd 'resonance' is actually a vibrational mode with a transfer function of the form ('s' being the Laplace variable)

H(s) = (s^2 + 2*zeta1*omega1*s + omega1^2) / (s^2 + 2*zeta2*omega2*s + omega2^2)

The denominator 'omega2' term corresponds to the peak at 4.99GHz, and the numerator 'omega1' term corresponds to the deep trough that follows the peak at 4.99 GHz.

One can demonstrate this behavior in Matlab very simply using the code below



zeta1 = .01; omega1 = 1.02; zeta2 = .011; omega2 = 1.0; 

num = [1 2*zeta1*omega1 omega1^2];

den = [1 2*zeta2*omega2 omega2^2];

freqs(num,den)



The frequency response is shown in the figure below which mimics your 2nd 'resonance'.


What you need to do for your 2nd 'resonance' is to estimate 4 (not 2) parameters, viz., omega2 and zeta2 (for the peak) and omega1 and zeta1 (for the trough). Depending on your data a nonlinear least squares optimization scheme should do. As a start, you could assume zeta1 and zeta2 to be equal (say zeta) and get an initial estimate of the terms 'omega1' and 'omega2' from your graph (or data). This leaves you with solving for just 'zeta'. Then you could iterate the model to get zeta1 and zeta2, etc., etc.

Note that the 'Q' term that you used is directly related to the quantity 'zeta' in my exposition here.