How can I estimate transfer function from frequency data?

Hello, I have magnitude and phase data of a system in frequency domain:

I want to derive its transfer fuction with the following code:

load('Magnitude.mat');

load('Phase.mat');

load('freq.mat');

data = frd(Magnitude.*exp(1j*Phase*pi/180),freq); % Frequency-response data model

np = 6; % # of poles

nz = 6; % # of zeros

iodelay = NaN;

sys = tfest(data,np,nz,iodelay);

bodeplot(sys,{2*pi*1e0,2*pi*1e6})

sys.Report.Fit

h1 = gcr;

setoptions(h1,'FreqUnits','Hz')

setoptions(h1,'MagUnits','abs')

Here is the result:

The estimated TF does not fit well with the phase angle but fit percent report is 98% and MSE is very high (3.6e6). I tried it for any number of zeros and poles but I got even worse results.

Could you please help me to solve it?

Thanks

[Mathieu NOE]

hello

it would be nice to share the data as well

[Easa AliAbbasi]

Hello @Mathieu NOE, thank you for your response.

I attached my data as well.

Answer 1

[Star Strider]

The best way to see how well the estimated system fits the data is to use the compare function.

Change (increase) ‘np’ using only the number of poles (initially use the default value for the number of zeros and if necessary, fine-tune later with ‘nz’) until you get an appropriate fit.

EDIT — (13 May 2022 at 14:50)

Increase the system order. I was able to get an acceptable fit to both the magnitude and phase with:

ord = 22;

sys = ssest(data,ord);

Also:

figure

compare(data,sys)

% sys.Report.Fit

h1 = gcr;

setoptions(h1,'FreqUnits','Hz')

setoptions(h1,'MagUnits','abs')

.

[Easa AliAbbasi]

Thank you. I increased the order to 12 which is the best order calculated by MATLAB for my system

and here is the report:

FitPercent: 99.0501

LossFcn: 1.2070e+06

MSE: 1.2070e+06

FPE: 1.9693e+06

AIC: 1.7322e+03

AICc: 1.7482e+03

nAIC: 14.4837

BIC: 1.7947e+03

The values of MSE and FPE are concerning me. Do you have any comments on it?

[Star Strider]

No.

Keep experimenting until you get the desired result.