How to estimate transfer function with findfit

Question

I already get a point series of the bode plot,(that is, system gain in frequency domain),{{freq1, gain1},{freq2,gain2}....} named as $data$

So I first assume that the transfer function is in the form of $$model=\frac{a1s^2+b1s+c1}{a2s^2+b2s+c2}$$ (by analyzing the bode plot, the tf should be in this form)

and then I use $FindFit$ function to estimate parameters $vars=\{a1,b1,c1,a2, b2,c2\}$

$$FindFit[data,Abs@model,vars,s]$$ however what I get is totally wrong, the fitting result gives a very high peak around 0 and never match the shape of original bode plot.

so, can you guys give some advice on how to estimate transfer function with mathematica

I know matlab has tfest. however, in my case, for some reason, it is better to do that is mathematica

\

Here is an example showing my process, with some simulation data, not the real experimental one.

fitmodel=Abs[n1/(a*(s*I)^4 + b*(s*I)^3 + c*(s*I)^2 + d*(s*I) + e)];
par={n1, a, b, c, d, e};
givenpar = Thread[par -> {500000, 1, 30, 26300, 270000, 25000000}];
givengain = Table[fitmodel /. givenpar, {s, 1, 300, 1}];
fitpar = FindFit[givengain, fitmodel, par, s];
fitgain = Table[fitmodel /. fitpar, {s, 1, 300, 1}];

Answer 1

The issue you have is with initial estimates (as noted by bill s). First I look at your data.

fitmodel = Abs[n1/(a*(s*I)^4 + b*(s*I)^3 + c*(s*I)^2 + d*(s*I) + e)];
par = {n1, a, b, c, d, e};
givenpar = Thread[par -> {500000, 1, 30, 26300, 270000, 25000000}];
givengain = Table[fitmodel /. givenpar, {s, 1, 300, 1}] // N; data = 
 ListPlot[givengain, PlotRange -> All]

This looks like two resonant modes with one resonance at about 25 and the other at about 150. To help getting good starting estimates I suggest you formulate your model in the form

model = a/Abs[( s^2 + b1 s + c1) ( s^2 + b2 s + c2)] /.s -> I ω

(* a/Abs[(c1 + I b1 ω - ω^2) (c2 + I b2 ω - ω^2)]  *)

Note that I don't have a coefficient on the highest power of s in the denominator. The OP does but this cancels with a factor of the term in the numerator and is not necessary. It is now easier to suggest starting values; c1 and c2 are the resonance frequencies squared we have already estimated. The values for b1 and b2 can be estimated from the sharpness of the peaks -I will take a guess. Now I fit using NonlinearModelFit which I prefer to FindFit because it gives good properties of the fit.

fit = NonlinearModelFit[givengain, 
   model, {a, {b1, 0.1}, {c1, 25^2}, {b2, 0.1}, {c2,150^2}}, ω];
fit["ParameterConfidenceIntervalTable"]

You can see that the parameters have been found with only numerical error. I am concerned that b1 is negative but this is the consequence of fitting the absolute values. If you have complex values you should fit those. Otherwise you may need to change the sign on the b values.

We can compare the model and the fit

Show[
 data,
 Plot[res[ω], {ω, 1, 300}, PlotRange -> All]
 ]

The model goes exactly through the data.

It is straightforward to convert my coefficients to yours but this is not the issue here. The problem is how to find good initial estimates. I often use some peak identification algorithm to find the frequency of the resonances and use these as the estimates.

Hope that helps.