I applied the prony with this now:
% % Prony method
[b,a] = prony(yt,16,16);
[r,p,k] = residuez(b, a);
s = log(p))*Fs;
t16 = t(1:16);
amp = r;
st = s'.*t16;
ht = amp'.*exp(st);
figure(4);
plot(t16, ht); grid on;

But this will give me h(t) of only 16 points and not matching to my response function. Now what should i do next.

 

I applied the prony with this now:

% % Prony method
[b,a] = prony(yt,16,16);
[r,p,k] = residuez(b, a);
s = log(p(:))*Fs;
t16 = t(1:16);
amp = r;
st = s'.*t16;
ht = amp'.*exp(st);
figure(4);
plot(t16, ht); grid on;

But this will give me h(t) of only 16 points and not matching to my response function. Now what should i do next.

It seems to me that you still do not understand the prony method. It looks like you tried the prony approximation with an order of 16. I have no idea if this is a good guess or not. However, assuming this is correct, why do you create a time vector with 16 points and then multiply each time value by one of the 16 poles. Is there some theoretical basis for doing this?

Instead, I would think you want to generate a larger time vector with an appropriate sample time. Then generate the approximate h(t) function as a sum of exponential terms with the amplitudes (amp) and frequencies (s) you calculated. I'm attaching the file I posted previously to remind you of the form of the Prony approximation. Here N=16, beta_i=amp(i) and s_i=s(i).

 

I analyze my system and found that it is of order 2. When i apply prony and residuez i got s and amp vectors of 3 points but my data and time is 500 points now how i multiply this and sum i.e. h(t) = sum(amp*exp(s*t)).
After getting h(t) what should i do toget my input from output, either u(t) = deconv(y(t), h(t)) or first take fft of h(t) and y(t) then u(f) = y(f)/h(f) and after this u(t) = ifft(u(f)).

 

Sir! steveb please give some response

 

Sir! steveb please give some response

I'm not sure why you get a vecor of length 3 when you have an order of 2. I think you have a mistake somewhere in the code. If you post your code, I can look.

Anyway, once you have the amplitude and frequency for each component of the summation approximation, you can take the inverse Laplace transform of the h(t) using a transorm table to get H(s). The inverse transform of an exponential function is well-known and available in all transform tables. Now you have the continuous time transform in equation form.

Now, the methods to get the input function are available. The ones you mention are possible and may work. However, you can probably just take the inverse of H(s) (i.e. 1/H(s)) to get the transfer function of output to input. Then lsim or other matlab tools should get you what you want.

 

Sir! my code portion:
%time t 500 points, ut input 500 points, yt output 500 points
%yf = fft(yt); uf = fft(uf);
h = yf./uf; % calculated t.f.

htt = abs( ifft(h,resolution) ); % calculated t.f. h(t)
order = 2;
[b,a] = prony(ht, order, order); % take prony of h(t) which %gives a and b vectors of size 3
[r,p,k] = residuez(b, a); % calculate r, p, k from a and b
s = log(p))*Fs; % si column vector
amp = r; % amplitude column vector
% Here i have problem in the next lines how it will be %calculated
st = s * t; % si * t
htc = amp.* exp(st);
ht = sum(htc); % h(t) calculated
hf = ifft(ht, resolution); % H(f) continuous form
% yt2 another output provided
yf2 = fft(yt2, resolution);%output in freq for yt2 (output)
uf2 = yf2./hf; % calculated uf2 (input) from yf2 and H(f)
ut2 = ifft(uf2, resolution); % reconstructed input from yt2

 

When I run your code, both "s" and "amp" are vectors of length 2.

Double check and try to find out why you get length of 3. My guess is that you first ran the program with order=3, and then re-executed it with order=2. If you do this without executing th command "clear" then the vector will stay at length 3 and will not be reset to 2.

It's sometimes a good idea to put the command "clear" at the top of your program.

 

Yes sir you are right that s and amp vector are of length 2. Now, i applied after getting the s and amp.

syms a1 a2 s1 s2 t w %symbols
mycont_tf = lapalce(a1*exp(s1*t) + a2*exp(s2*t)) % gives laplace
% mycont_tf = a1 / (s - s1) + a2 / (s - s2)
Since i have amp and s vectors so i get mycont_tf. Is it correct. Now, if i apply any of yt i can get ut.

 

Yes sir you are right that s and amp vector are of length 2. Now, i applied after getting the s and amp.

syms a1 a2 s1 s2 t w %symbols
mycont_tf = lapalce(a1*exp(s1*t) + a2*exp(s2*t)) % gives laplace
% mycont_tf = a1 / (s - s1) + a2 / (s - s2)
Since i have amp and s vectors so i get mycont_tf. Is it correct. Now, if i apply any of yt i can get ut.

Im not sure. I think that you now have the transfer function H(s) for input to the output, but to get the TF from the output to the input will be the inverse function 1/H(s). I can't be sure because you may have solved for the inverse function directly.

After that, I dont' know if it's correct for sure.

Proving out any computer program is not easy. When I wrote my code for solving my problem, I first applied the program to test situations where I knew the answers ahead of time. Basically I would make a simple system in Matlab (second order for example). I would calculate the system poles, transfer functions etc. Then I would apply input signals and calculate output signals. Then I would apply my program to the data and see if I could extract the impulse response correctly.

This is the type of thing you need to do. Run test cases under ideal conditions (no noise) with systems that you know exactly. Then see if the program is working. If it works in the ideal cases, then add some gaussian noise to your data and apply the method again to see how much noise can be present before the method fails.

Once you go to the real data, you will have to use judgement and decide if the noise in your data is within acceptable bounds. If the noise is too large, you will need to go to more sophisticated methods.

Basically I can't do your work for you, but these basic guidelines seem reasonable to me. As I mentioned early on, this is not an easy problem.

 

I have a new idea for one of my old job. As you know that: At the input i have a signal x(t) from an device sampled at 10 nSec at scope (500 data points) and that signal is communicated through a coaxial cable (approx 1 Km length) which due to attenuation attenuate that signal and recieved a sampled attenuated signal y(t) (sampled at the same rate as at input ). I have only one input signal/vector but online ouput updating. My objective is to reconstruct the original signal at input from the signal recieved at output. This system is online (continuously updating signal at ouput) but i have only one input. For this i guess that cable will be my system and i need the impulse response h(t) of it so that i deconvolve my output y(t) with it to get my input x(t). For time-domain deconvolution i want to use conjugate-gradient method (Deconvolution of Impulse Response from Time-Limited Input and Output: Theory and Experiment ) first to get impulse response h(t) from input x(t) and output y(t) deconvolution and then use it h(t) to get input x(t) from impulse response h(t) and output y(t) deconvolution. Am i right or not. I have no idea of how to implement it. Please! comment and help.