I have input data u(t) and output data y(t) of a system having 500 points. This is test data from a SISO system. I am trying to use MATLAB so that it take a vector of input data, output data, and a sample time and produce the frequency response of the system (magnitude and phase). After this i want to make transfer function for this. My objective is to reconstruct the input signal from my output or output from input. Actually i need reconstruction signal of my system and for that i think i will use that transfer function. i.e. in lsim.
I use a manual method of finding poles and zeros of the response. How, i can calculate the zeros and poles from the response data. Am i right to reconstruct the signal input from output and output from input. Is there any other way to reconstruct such signal. How can i get T.F. from my calculated response. If anyone understand and can solve my problem please! help. How can I do this in matlab or other method.

 

Do you have access to the following Matlab toolbox?

http://www.mathworks.com/products/sysid/

 

I have input data u(t) and output data y(t) of a system having 500 points. This is test data from a SISO system. I am trying to use MATLAB so that it take a vector of input data, output data, and a sample time and produce the ..... transfer function for this. ....... How can i get T.F. from my calculated response. If anyone understand and can solve my problem please! help. How can I do this in matlab or other method.

Reconstruction of the transfer function from the input and output samples is generally not a trivial thing to do, but there are methods. I've done this using the Prony method which relies on fitting the transfer function to a sum of exponential functions. There are many issues in doing this. If the data is taken from experiments, then noise can prevent the simple method from working and then you need more sophisticated techniques. If the data comes from calculations and is noise free, the simple Prony method works well, but is usually approximate, but usually you can identify the correct number of poles and estimate the values.

I would only recommend going down this path for serious work like a college semester-long project, graduate thesis or a company work project.

The following link would just be a starting point. If you do go down this path, I can dig up some references I have at work. I might even have an old Matlab file which illustrates the process.

http://en.wikipedia.org/wiki/Prony's_method

Matlab has two functions (probably in one of the toolboxes) useful to do this task; prony and residuez. The "prony" function calculates the Prony coefficients and the "residuez" function does a partial fraction calculation to help extract the poles.

 

Tell me more about prony method. May be i got success. I need reconstruction of the transfer function from the input and output samples.
Uptill now i have done this in MATLAB:

clc;
clear all;
t =0.4e-9*(1:1:500);
ut = [500 points];
yt = [500 points];

resolution = 4096;
a1 = resolution/2;
Fs = 2.5e9;
uf = fft(ut, resolution);
yf = fft(yt, resolution);

h1 = yf./uf;
h = h1(1:a1);
f = Fs * (1:a1)/resolution;
gain = 20 * log10(abs(h));

figure(1);
semilogx(f, gain); grid on; hold on;
axis([1e5 1e8 -30 0]);

p = 1;
% __________________________________
n1=conv([p/5.25e7 1],[p/5.3e8 1]);
d1=conv([p/3.35e7 1], [p/9.5e7 1]);
%___________________________________

sys1=tf(0.72*n1,d1);
out1=lsim(sys1,ut,t);

n2 = d1;
d2 = n1;
sys2=tf((1/0.76)*n2,d2);
out2=lsim(sys2,yt,t);

bodemag(sys1); grid on;
bodemag(sys2); grid on; hold off;

figure(2);
plot(t,ut,t,yt,t,out1,'r'); grid on;

figure(3);
plot(t,ut,t,yt,t,out2,'r'); grid on;

In this i have got pole & zeros manually by response plot. Is there any proper method so that i can get these and reconstruct my signal. Sysid toolbox can't help me in this regard becuase those poles values doesn't match to my values and i have no other idea about it. Can anyone help.

 

Tell me more about prony method. May be i got success. I need reconstruction of the transfer function from the input and output samples.
Uptill now i have done this in MATLAB:

I can give some general information based on what I did in the past, but keep in mind that you will probably have to do some of your own digging too. Your problem may have some particular issues, such as signal noise, that require further considerations.

First, here is another web-site that might help:

http://www.statsci.org/other/prony.html

Second, I made a PDF for a simple overview of the method.(prony_method.pdf)

Third, I've attached a paper that summarizes the more general problem and talks about issues and solutions if noise is present. This general problem is difficult and an ongoing active area of research. (noise.pdf)

Forth, I extracted the critical lines of code from my Matlab file as follows:

order = 16;                        % order of prony approximation
[b,a]=prony(DATAFILE,order,order); % calculate Prony coefficients
[r,p,k]=residuez(b,a);             % do partial fraction expansion
s=log(p(:))/Ts;                    % extract poles in real frequency
amp=r;                             % extract amplitudes 

This code looks very simple, and it shows how powerful some of the Matlab tools are. But, it took me a couple of days to understand the method and make the connection that these tools would work. If you want to implement something like this, understand the method, then understand those Matlab functions.

 

Thank you i will try it and trhen respond. You didn't tell about reconstruction. What about that.

 

You didn't tell about reconstruction. What about that.

Do you mean reconstruction of the transfer function in the Laplace frequency s-domain? If so, the method above allows the system impulse response to be approximated by a sum of exponential functions. You can then take the Laplace transform to get the transfer function in the frequency domain. Laplace transformation of a sum of exponentials is straightforward using a Laplace Transform table. This allows extraction of poles and zeros. So basically the system transfer function is reconstructed in equation form from the time-sampled input and output signals.

The question as to whether you can accurately determine the correct number of poles and zeros and provide reasonably accurate estimation of the values of the poles and zeros will depend on how clean the time-domain samples are, and how good your intuition is in applying the method. You must guess the number of exponentials to use and then see how good the fit is. The quality of the fit of the approximation to the real data will make it obvious when the number of terms is too small. Generally you want to use just enough terms to get a good estimation.

Keep in mind that I've provided just an outline of the method. Many things will become obvious if you start to study this approach and implement it for your specific problem.

 

I want to hear about system signal reconstruction. Am i right in the method, if not then how can i reconstruct output from input or input from output (i have time domain 500 points data of input and output). Please tell me the steps.

 

Sysid didn't help me.

 

I have applied prony and residue. I also use stmcb ( insted of prony) because that is better than prony. Now, i have a question that through residue i got P.F. how can i get polynomial (numerator) from that because i need numerator and denominator polynomials. Poles give me the denominator polynomial ( through Poly(p) ). Give some idea about it.

 

I have applied prony and residue. I also use stmcb ( insted of prony) because that is better than prony. Now, i have a question that through residue i got P.F. how can i get polynomial (numerator) from that because i need numerator and denominator polynomials. Poles give me the denominator polynomial ( through Poly(p) ). Give some idea about it.

I'm a little confused as to what you are asking, but I'll take my best guess.

I wasn't aware of the "stmcb" function, but that looks very useful, and I will look into this next time I have this type of problem.

Both "stmcb" and "prony" return the numerator and denominator coefficients of a digital system that matches the sampled impulse response (or sampled output/input data).

If you only need the digital system, then you are done, but I was assuming that you were interested in a continuous time transfer function that would have produced your sampled data. However, looking back I see you never said whether your system was digital or continuous.

In the method I recommended, using "residue" and the data extraction, you get the continuous time transfer function as a partial fraction expansion (this is a sum of exponentials in the time domain). At that point, you can just add up all the terms using a common denominator and then multiply out the numerator and denominator to get the transfer function as a ratio of polynomials (in terms of complex Laplace frequency "s").

Did I understand your question? If not, please tell me if you are trying to reconstruct a digital system or a continuous system from your data.

 

Yes i want continous system transfer function from my first time input-output data so that when next time i provide any output to it i get the input.
I have some problem in getting numerator from P.F.

 

Yes i want continous system transfer function from my first time input-output data so that when next time i provide any output to it i get the input.

OK, good. The method I showed enables you to do this.

I have some problem in getting numerator from P.F.

I'm probably being absent minded, but what does P.F. stand for?

The code section I posted previously shows how to extract the continuous time complex amplitudes and frequencies for the Prony approximation, which is just a finite summation of exponentials. This sum of exponentials in the time domain can be transformed into the frequency domain and the result is the partial fraction expansion of the transfer function. In general, if you have a partial fracton expansion you can combine all the terms into an equivalent transfer function which is a ratio of polynomials. You just find the common denominator and do the math.

If this doesn't make sense, then I suggest you post your work and I can try to figure out where you are in the process.

 

I mean from P.F (Partial Fraction). I have attached my work done in MATALB. Guide me for this in signal reconstruction.

 

I mean from P.F (Partial Fraction). I have attached my work done in MATALB. Guide me for this in signal reconstruction.

 

IGuide me for this in signal reconstruction.

Honestly, I can't understand what you are asking. Before you said you were having trouble getting the numerator from the PF expansion. Now you are saying to guide you on the signal reconstruction.

The matlab file you gave me is not commented well and I don't have time to figure it out. When I run the program I get three graphs which I don't understand, and they don't seem to indicate that you have found a transfer function that matches your data.

At this point I'm not sure exacly what you are trying to do, or where you are in the process. I really can't help much further other than to give the basic approach as follows.

1. Create a data vector which has the samples of the impulse response function of your system. (Do you know what an impulse response is? The input signal convolved with the impulse response function gives your output signal)

2. Use the "prony" function with the data vector and a guess at the order of the system .... i.e. use the following:
prony(DATA,order,order). This will give you the numerator and denominator coefficients for the digital system transfer function.

3. In order to get to a continuous system, do a partial fraction expansion using "residuez". This form is useful because if you transform each term back into the time domain you will get a sum of exponential terms that can be matched up with the Prony approximation as a sum of exponentials. (Do you understand the Prony method? You must study and understand this method thoroughly for this step to make sense)

4. Use the code I gave you to extract the amplitudes and frequencies of the exponential terms and convert to an equivalent continous system representation, as shown below. This is a bit confusing and you may just want to take my word that this works. Note that you need your sample time for the data which is shown as Ts below.

DATA=????? % create vector of your impulse response samples
order=4; % make a guess at the order
[b,a]=prony(DATA,order,order); % calculate Prony coefficients 
[r,p,k]=residuez(b,a); % do partial fraction expansion
s=log(p(:))/Ts; % extract poles in real frequency
amp=r; % extract amplitudes

Remember this is an approximate method. How well it works depends on whether you have low-noise data and if your system really is a low order system that can be represented by the form of the Prony approximation. Even then, you have to guess the order and possibly limit your data set if it decays near zero at the end.

There is a bit of an art in doing this. As I said, it took me a few days to get this to work in my application. Before you try you real application, try a test case with a known simple low-order system. For example, make a transfer function of a second order system and use matlab to create a sampled version of the impulse response in the time domain. Then you can try to reconstruct the original transfer function from the sampled data.

 

In fact i have no impulse response function, only have input signal and output signal data (500 points each) with a time vector.
From these vetors i need an equation so that when i provide another output signal, that give me input signal. The work done so far is attached in file with comments (actually previous file has some missing data points).

 

In fact i have no impulse response function, only have input signal and output signal data (500 points each) with a time vector.
From these vetors i need an equation so that when i provide another output signal, that give me input signal. The work done so far is attached in file with comments (actually previous file has some missing data points).

OK, I'll look at this new file. I'm busy with a few things at work, but maybe tonight I'll look at it.

In the mean time. I understand what you are trying to do and the method I showed is one way to get there. You can construct a sampled version of the impulse function directly from the input and output data. With that you can use the Prony method to get an estimated continuous time transfer function. Once you have that, you can generate input data from any other output data.

There are a several different ways to do these various steps, either in the time domain or in the frequency domain. I'll try to think about the best method for your situation, but some of these steps are commonly done and perhaps some other people can provide suggestions. For example, the following two steps have several possible methods. What are the most accurate and simplest to implement?

1. Generate a sampled impulse response function h(t) from the sampled input and output data.

2. Given a known continuous time system, generate input data from output data.

 

1. Generate a sampled impulse response function h(t) from the sampled input and output data.

2. Given a known continuous time system, generate input data from output data.

I didn't have a chance to put much time into this. I noticed that you used the fft to generate the transfer function. You could use this and transform back to the discrete time domain. You would then have a sampled version of the impulse response function (approximately).

You can then use the Prony method to get the continuous time version of an approximate system transfer function. Then invert this to get the inverse system. You can then use "lsim" with this inverted transfer function to get the input signal from the output signal.

I'm just throwing ideas out right now. I'm not sure any of my suggestions are the best method, but I think the approach should work.

 

The work done so far is attached in file with comments (actually previous file has some missing data points).

In looking at your code I see that you have misinterpreted something about the code that I gave you before. Referring to the section of code shown below. Why did you chose to make the denominator of your transfer function (den) equal to the poles? Also, why are you even trying to find a numerator and denominator at that point? This shows that you do not understand the Prony method. My original example told you to extract the poles and the amplitudes needed in the Prony approximation. The Prony approximation I provided in one of my first posts. It is a time domain approximation of the impulse response as a sum of exponential functions. Once you have the complex frequencies and the complex amplitudes of the exponential terms, you have your system impulse function in continuous time.

% [b,a] = prony(h1,6,6);
% [r,p,k] = residue(b, a);
% den = log(p(:))*Fs;
% %num = ????? Here is i confused.
% sys = tf(k*num, den);
% input = lsim(sys, yt, t);