How to estimate transfer function from real measurements data (Partially known)

[Nemesis]

I will put it here again.

[WCR]

Nemesis:

You may want to look at how I attacked a similar problem, as published in Am J Physiol 265: H1577-H1587, 1993.  I had two inputs (cardiac inflow & outflow) and two outputs (aortic and vena cava pressure), so I had to estimate 4 transfer functions (4 impedance spectra).  You have a single input and single output, therefore just one transfer function to estimate.  First I estimated the 4 transfer function spectra in the frequency domain using standard statistical techniques. I used a "nonparametric" approach, i.e. I used FFTs, rather than making a "parametric" (ARMA-type) model.  Then I identified the RLC circuit model (which is what you said you would like to do) which gave the best fit to the 4 spectra.  "Best fit" to the 4 spectra was defined in a statistical way.

[WCR]

Like Thierry C, I cannot open any of the attachments to your ealrier posts.  I can open your most recently posted files: results.txt and Signal_construction.vi.  What are columns 1, 2, and 3 in results.txt?  Column 1 = frequency in radians/s?  (I notice that col. 1 has equally spaced values from 1*(2Pi/10) to 100*(2Pi/10).)  Col 2=gain in dB?  Col 3 = phase?  Or maybe something completely different.  And if results.txt IS a transfer function, what is it the transfer function of?  Pressure divided by flow?  Pressure divided by velocity?  Pressure divided by piston displacement?  Other?

Bill R.

[RavensFan]

There was probably nothing wrong with those attachments to begin with.

There has been some problems with the forums lately killing some attachments and other things such as signatures.

Problem with Attachments

Hopefully once the problem is debugged, the attachments can be restored.

[WCR]

You said "I believe and I know that the system is hardly linear. In my system I have pipes that are not rigid. They are soft such as a silicone pipe. If we simple compare hydraulic system to the electric one we can write that it is like a RLC system for each of a pipe (a picture).  We have friction effect, viscous effect, compliance effect and viscoelastic effect of pipes. Then 7 parts of these kind of system with different resonance plus another equipment sych as pumps, valves etc."

A hydraulic system with elastic or viscoelastic pipes can be treated as linear, for a range of flow regimes.  In other words, the Navier-Stokes equations assume a linearized form, if certain assumptions hold.  This was demonstrated by John Womersley in the 1950s.  You'll want to become familiar with his work if you are not.  Some of his best papers were techinical reports for the US Air Force and are a bit hard to get and hard to read, so you might want to read the summaries of Womersley's work in "McDonald's Blood Flow in Arteries", by Wilmer Nichols and Michael ORourke, or see the book "Hemodynamics" by the late and esteemed Bill Milnor.  Or see these papers by me and others:

Johnson, D.A., Rose, W.C., Edwards, J.W., Naik, U.P., Beris, A.N., Application of1D blood flow models of the human arterial network to differential pressure predictions. J. Biomech. 44: 869-876, 2011.

Johnson, D.A., Spaeth, J.R., Rose, W.C., Edwards, J.W.., Naik, U.P., Beris, A.N., An impedance model for blood flow in the human arterial system. Part I: Model development and MATLAB implementation. Computers and Chemical Engineering, 35: 1304-1316, 2011.

[Nemesis]

Thanks for answer WCR. It is nice to hear that somebody shares similar work. You wrote three post I will answer to every of them here below.

Post 1

As you said at that moment I am trying to get a good flow signal. After that I want to combine both of in vivo signals i.e.. flow and pressure. With my previous system I just looked at the difference between measured gain or/and phase and estimated ones and I checked how much it changed per decade. 20dB or 40dB etc. Then I knew it if I have to add first or second order function to my existing transfer function. Recently I changed the pipes and length and the function change as well. Unfortunately for me, I cannot get a good transfer function with my previous technics. Therefore I am searching help in Labview functions. Just one thing you consider your system as a single RLC system (Lumped model) or something more complicated? Normally for me it is several RLC systems. Anyway, I will read your article.

Post 2

In the results.txt. In columns:

1 - Pulsation in [rad/s]

2 - Gain in [dB]

3 - Phase in [deg]

The length of measured frequencies is 0.1 Hz with 0.1 Hz step to 10Hz. And it is the transfer function for it:

num = 1.226e-006 s^4 + 0.0001119 s^3 + 0.00206 s^2 + 0.07692 s + 0

den = 5.676e-011 s^7 + 3.633e-009 s^6 + 4.23e-007 s^5 + 1.453e-005 s^4 + 0.0006311 s^3 + 0.01016 s^2+ 0.1717 s + 1.

The transfer function is the movement of the piston in V to flow measured on the flowmeter in V. Or as you maybe prefer flow divided by piston displacement.

I added new data in the attachment that I have problem and I cannot find estimation.

Post 3

I know those people that you mentioned about. My flow calculation are based on Womerslay's equations after I check with CFD calculation as well. I based on his and their articles. I have not read "McDonald's Blood Flow in Arteries" yet. Generally I based on quite good books such as: 

M. Zamir The physics of pulsatile flow.2000.

Y.C Fung . Biomechanics. Circulation. Second edition 1996

Actually, I met some opinions that I should go further and not to look at the something made in 1950 ... but I think it is still a great tool. I think, we could write many posts speaking about our work but I think it is not the case of this forum. Anyway I am open to any discussion. Once again I will check your articles.

Kamil

[Nemesis]

Once again. Basically you wanted to put all in the similar equations (look at the attachment). However it is flow to pressure or opposite. I need to have motion of piston to flow. I can assume that speed of motion of the piston dh/dt multiply a surface of the piston give me flow (volume in time) that I add to my system. Eventually I can assume that force with surface give me pressure. Hmmmmm I have to take a closer look at this but I do not know if I can get something sufficient for estimation.

Kamil 

[WCR]

I understand better now. Your transfer function is flow q(t) (measured by flow probe) divided by piston displacement x(t).  The transfer function you have measured reflects the properties of the displacment measrurement (how well the "displacement voltage" represents true piston position), the flow measurement (how well the "flow voltage" reflects true flow) and the hydraulics between the piston and flow sensor.  If the piston displacement and flow signals are accurate, and if the piping between the piston and flowmeter is rigid, and if the water is incompressible and does not cavitate, then the flow should be the time derivative of displacement: q(t) = A*dx(t)/dt, where A=area of the piston.  Which you knew already.  Therefore we would expect a transfer function G(s) = Q(s)/X(s) = A*s, which would have a slope of +20 dB/decade and a phase angle of +90 deg at all frequencies. 

The transfer function you measured (results.txt) resembles the ideal transfer function at some frequencies, although it is not perfect, of course.  You have estimated the coefficients for a rational polynomial approximation to your measured transfer function:

G(s) = (b1*s + b2*s^2 + b3*S^3 + b4*s^4) / (1 + a1*s + ... + a7*s^7)

You noticed that when you changed the tubing attached to your system, the transfer function G(s) changed.  This may or may not be a surprising result, depending on the particulars of your system.  

You would like to know if there is an RLC circuit (or hydraulic analog) that has a transfer function similar to the one you measured, and if so, what are the values of the circuit elements.  You will have to calculate the theoretical transfer functions of different circuits with different numbers of Rs, Ls, and Cs, and adjust R,L,C to get a good fit to your measured transfer function.  Labview, Matlab, etc., have tools to help with this multidimensional minimization. I recommend fitting the measured transfer function (the data in results.txt), rather than fitting the rational polynomial approximation, which is itself a fit of the measured data.  Of course your "goodness of fit" function, or "error function", should somehow take into account the fact that you may have more confidence in your transfer function as some frequencies than at ohthers.  For example, it is common to have larger "error bars" on the transfer function estimate at the highest frequencies.  Also, you need to considert how to do the dfitting of complex data: match the magnitude, or the phase, or both?  I recommend minmizing distance in tthe complex plane, i.e minimize the sum of the error in the real part and the error in the imaginary part, added up across all frequencies.

I agree with you about not looking at papers that are from 1950!  Fung and Zamir are excellent sources.

[Nemesis]

Hello,

Thanks for answer. I did addition measurements today. I put in attachment. I made experiments for different of amplitudes. As you can notice there are some differences in the range of frequency 0.1-15Hz. Therefore I have some difficulties to get nice signal (with this configuration) and I think I need quite accurate transfer function. I will fight still. I will take your directions into consideration.

Kamil

[WCR]

If it were me, I would try including a pure delay in the transfer function (a factor of exp(-s*T) in numerator).  The linear nature of the phase response at high frequencies suggestes that a delay is present.  The magnitude of the delay, T, can be estimated from the slope of the phase curve at high frequencies:

T~=(-slope of phase response, in degrees per Hertz)/360.