Using ARX to obtain a transfer function from frequency response data

Hi,

I have some complex frequency response data obtained from a measurement. I'm trying to fit a transfer function to this data using the arx function from the system identification toolbox.

However, when I compare the calculated model to the measured data, they don't match very well. Using some knowledge about my system I can calculate the coefficients by hand and this transfer function matches the measured data perfectly. So I'm sure that the system can be approximated by a transfer function of the specified order. As the order of the system is not very high and the measured data is of good quality I think the identification should be better than the results I get, but I can't figure out what I'm doing wrong.

I have included the commands I'm using below. Any help would be very much appreciated. Thanks in advance, Christian

frdobj = idfrd(complexRespData,frequencies,0, 'Units','Hz', 'InputName','In', 'OutputName','Out', 'Name','Z', 'spectrumdata',[]); % create idfrd object from measurement data
modelobj = arx(frdobj,[2 1 0]); % identify. Use 2 poles and one zero.

modeltf = tf(modelobj);   % transform to tf
[ident_mag, ident_phase] = bode(modeltf,2*pi*frequencies); % obtain frequency response
model_comp = ident_mag .*exp(1j*ident_phase * pi / 180);
% compare the result to the measured data
figure(1);
loglog(frequencies,abs(complexRespData));
hold on;
loglog(frequencies,abs(model_comp));

Answer 1

[Rajiv Singh]

ARX is usually not a good estimator of frequency response data. Please try Output-Error (OE) model which more closely represents a transfer function: y = B/F u + e. Use the OE command for the corresponding estimation as in:

 modelobj = oe(frdobj,[1 2]);

You could also try state-space estimation using PEM, as in:

 modelobj = pem(frdobj, 2);

If you are using release R2012a of MATLAB, the recommended estimation commands for frequency response data are TFEST (transfer function model), PROCEST (low order transfer function expressed in pole-zero form), and SSEST (state-space model). See:

http://www.mathworks.com/help//ident/ref/tfest.html

http://www.mathworks.com/help//ident/ref/procest.html

http://www.mathworks.com/help//ident/ref/ssest.html

[Christian]

Thank you for helpful reply. I've tried both PROCEST and TFEST on my data.

While TFEST works considerable better than my arx approach, I've still got some trouble to get a good fit at low frequencies. I've tried to use the focus interval from tfestOptions to emphasize the mismatching frequencies, but to no avail.

At

http://postimage.org/image/s8w601nad/

you can see an image of the resulting fit. The blue line is the measured response, while the red dashed line is the tfest result. The green dashed line is the manually determined transfer function, which has the same order as the identified result.

The upper graph shows the magnitude, while the lower graph contains the phase.

I've also tried to use a weighting vector for focus, but it seems that the required length does not equal the size of the input data, so I'm a bit unsure how to use it.

I hope you can help me again to improve my results.

[Rajiv Singh]

Difficult to say more without looking at the actual data. The difference in fit looks more pronounced than it actually is owing to the log scale used for the magnitude (e.g., look at this in absolute/linear scale). You can check modelobj.Report.Termination to see why the estimation stopped. It it ran out of iterations, you could try increasing maxiter (opt = tfestOptions; opt.SearchOption.MaxIter = 100; modelobj = tfest(frdobj,2,opt), and also play with various search methods (e.g.: opt.SearchMethod = 'lm') )

[Christian]

Thank you again for your insight.

Changing the search method does not improve the result. You are right that the logarithmic scale makes the differences more pronounced, but unfortunately the relative error also matters for my application.

My approach has been to put some emphasis on the the lower frequencies by using

opt.Focus = [40*2*pi, 1e4*2*pi];

so that frequencies from 40Hz to 10kHz are used.

This improves the fit for this part dramatically but gives overall poorer results (0.2649%, FPE: 6.11538e-07 vs. 99.45%, FPE: 40.8088 for the full fit).

I'm now trying to use the focus weights vector but I can't get it to work.

I tried

opt.Focus = 1./abs(respData);

Perhaps you can tell me how to use it properly. No matter what length I supply, I always get an error that it does not match the input data length. I've tried the length of the frequency vector as well as 2*length(frequencies)-1 as I saw that complex conjugate is appended to the data during the estimation. Both lengths result in an error at different stages of the estimation.

[Rajiv Singh]

If you have single-input data, try this:

fd = complex(iddata(frdobj))

Now prescribe weights for output samples (fd.OutputData) as a column vector of length equal to lenght(fd.Frequency). Use fd as estimation data.

This is a not convenient workflow for weight assignment and will be improved in forthcoming release.

[Christian]

Many thanks for your patient and helpful replies. Thanks to you I'm now able to get a good fit even at lower frequencies. I'm documenting the method I've used below, just in case someone else might have similar demands in the future.

The way it worked for me was to estimate a model, specifying only the number of poles and zeros and leaving the default options.

modelobj1 = tfest(frdobj,N_poles, N_zeros,0);

I then use your method to give all data points the same weight, regardless of their absolute value by using

opt = tfestOptions();

fd = complex(iddata(frdobj));

weights = (1./abs(fd.OutputData));

opt.Focus = weights;

Then I use the first estimation as an initial guess and estimate a second time using the weights:

modelobj = tfest(fd, modelobj1 , opt);

The result is a transfer function with complex coefficients that shows a good fit on a logarithmic scale. It minimizes relative deviations instead of the absolute ones.

I still have to find out what meaning the complex coefficients have for my application, but the estimation is now stable for different samples of my input data.

Answer 2

[Li Zhijun]

I have been faced with the same problem as you. I followed your method with same weight, the quality of the fitting is greatly improved. But what do the complex coefficients mean is beyond my knowledge. I am looking forward to further discussion with you and i will search for a better solution or explanation.

[Christian]

I'm myself not quite sure if I fully understand the meaning of the complex coefficients. I'm not an expert in system theory, but as far as I can tell the complex coefficients effectively increase the system order. (try to make the coefficients in the denominator purely real by multiplying with the complex conjugate and you can see that it's in fact two systems in parallel: one with purely real coefficients, and one with purely imaginary coefficients in the numerator. The outputs of those two systems are summed. The system order is doubled.)

With my data, the imaginary coefficients are very small compared to the real ones for the same exponent. As my data comes from a complex impedance measurement I can attribute this imaginary system to 'second order effects', such as skin effect etc, which I didn't take into consideration when determining the number of poles and zeros for the identification.

If the real and imaginary parts of your coefficients have the same order of magnitude, I'd recommend that you try to increase the number of zeros and poles for your estimation. Hopefully, the higher you select your system order, the more the real parts will dominate, up to the point where you can neglect the imaginary part.

If this does not work you might want find someone who is more familiar with the details of system theory. Perhaps it's even possible to transform this transfer function into one with purely real coefficients analytically.

If you get a better explanation for the meaning of this complex coefficients, I'd be very interested to hear about it.

Answer 3

[lexi11]

Hi,

I am trying to follow your method but when I try the following, I get an error:

data = iddata(signal,[],Ts);
model1 = tfest(data,2,0);

error:

Error using tfest (line 107)
Index exceeds array bounds.
Error in test (line 14)
model1 = tfest(data,2,0);

My system is unknown. I am working on videos so my output is a signal that has pixel brightness values (one value for each frame for one colour channel). The signal is a single colour channel. And I am more interested in finding poles (I don't know how many poles or zeros it has). I tried with indicating zeros as 1 also, same error appears. I don't know what is wrong here.

Any help is appreciated.

Thanks.

[rupprechtt]

Hi,

did you solve the problem? I have the same problem.

Error using tfest (line 107)
Index exceeds the number of array elements (2).