How to find proper Impulse response function of the system with Gaussian signal input/output

[OpiCAN00]

I am trying to identify system by searching of proper impulse response function (IRF). I tested few methodology without succeed.

Source https://www.physicsforums.com/threads/transfer-function-laplace-of-the-impulse-response.783541/

From the basics of the signal theory the impulse response can be calculate like (L=Laplace transformation): L[x(t)]= X(s) , L[y(t)]= Y(s) and transfer function H(s)= Y(s)/X(s), then IRF is inverse Laplace transformation of transfer function  h(t)=L−1[H(s)]. By convolution of input signal (stimulus) and IRF h(t) can be calculate output (response) of the system y(t)=h(t)*x(t).

Since there is no inverse Laplace transformation function in LabVIEW and Laplace transformation function can handle just array of size smaller then 100, I tried to use FFT function instead. Discrete Laplace transformation use FFT, but it is not the same and probably there is my problem. I try to check, if impulse response is correct by convolution with input signal, but the result is never the output signal used to find IRF. I tried more ways how to find IRF but, as You can see in attached VI, non of it gives original response signal when convolve with input signal. I assume it is necessary to use some normalization of FFT, IRF and convolution also.

Impulse Response Function PtByPt VI and Deconvolution PtByPt VI of stimulus and response signal giving complete mess as a result.

SI Estimate Impulse Response VI give some result but after convolution with input signal, the amplitude of response is too high. The same result if I use FIR Filter VI  instead of convolution.

Direct descrete Transfer Function Estimation is too sophisticated and i am not sure about orders of numerator and denominator.

Please, if someone understand signal theory in LabVIEW tell me some way how to find IRF, which gives response signal after convolution with input signal.

Probably special solution is necessary, since In/Out signals are non periodical gaussians.

Thank You for answers..

P.S.: I tried a lot, to find solution in forums and study theory. 🙂 

Y(s)/X(s) = H(s)

Source https://www.physicsforums.com/threads/transfer-function-laplace-of-the-impulse-response.783541/

Y(s)/X(s) = H(s)

Source https://www.physicsforums.com/threads/transfer-function-laplace-of-the-impulse-response.783541/

Y(s)/X(s) = H(s)

Source https://www.physicsforums.com/threads/transfer-function-laplace-of-the-impulse-response.783541/

Y(s)/X(s) = H(s)

Source https://www.physicsforums.com/threads/transfer-function-laplace-of-the-impulse-response.783541/

Y(s)/X(s) = H(s)

Source https://www.physicsforums.com/threads/transfer-function-laplace-of-the-impulse-response.783541/

Source https://www.physicsforums.com/threads/transfer-function-laplace-of-the-impulse-response.783541/

X(s)=L{x(t)}(s)

Source https://www.physicsforums.com/threads/transfer-function-laplace-of-the-impulse-response.783541/

X(s)=L{x(t)}(s)

Source https://www.physicsforums.com/threads/transfer-function-laplace-of-the-impulse-response.783541/

[OpiCAN00]

I found partial solution for my problem. I took advantage of Gaussian symmetric property and did FFT only on half of the function starting from the center. I think the problem is that: If I do FFT on whole gaussian, the result include harmonics from gauss peak itself and impulse response function result in nonsense. Proper solution would be probably some kind of window or filter of FFT output but i am not relay sure how to do it. So my solution work just for symmetric functions and application of IRF for real input signal is more complicated but possible (divide signal in the center of the mass on halfs, then convolve one half regularly and second half solve rotated, result of convolution rotate back and concatenate results).

Still i would appreciate, if someone will share experience and tell me how to do it properly. 

Since  Gaussian is analytical function, its Fourier transformation can be found analytically. I did some calculations but discrepancy of result is too high. Possibly i did some mistakes in calculations or it have to be solved by Laplace, but Laplace transformation of Gaussian result in error functions and algebraic solution is very complicated.

stimulus:
f(t)=(A/(s Sqrt[2 Pi])) Exp[(-t-c^2)/(2 s^2)]

response:
g(t)=(B/(r Sqrt[2 Pi])) Exp[(-t-d^2)/(2 r^2)] 

where:

A- Aera of stimulus
B- Aera of response
c- Center of stimulus
d- Center of response
s- Sigma of stimulus
r- Sigma of response 

than impulse response is:

h(t)=F-1{F[g(t)]/F[f(t)])}= B/A Sqrt[2]/(E^((t+c-2d)^2/(4 (r^2 + s^2))) Pi Sqrt[r^2 + s^2])

[OpiCAN00]

I tried another method how to find IRF. Thanks to knowledge that IRF of system with Gussian Out/In has to be also Gaussian, is possible to find IRF by fitting response by formula including stimulus convoluted with IRF. Result of fit is parameters of IRF gaussian (fitting just require good guess of initial parameters of IRF). I can Reconstruct response by convolution of stimulus and IRF with high precision. Problem appear, when i try to reconstruct signal from two stimulus peaks in one spectra. Result of convolution have shifted second peak and i do not know how to normalize calculations to receive result from convolution with non shifted peaks. I figured out that linear convolution output array size is equal to Xsize+Ysize-1, based on that knowledge i did some normalization of X axis. But I still struggle with that task a lot. I expected some behavior of convolution result like You can see on the attached image.  

Convolution

But i receive this shifted peaks

Please can someone explain how to handle this?

[OpiCAN00]

Finally i figured out how to normalize convolutions outputs. Thanks to this post : https://forums.ni.com/t5/LabVIEW-Idea-Exchange/Output-size-configuration-for-1D-convolution/idi-p/93...

Polymorphic version of convolution of 2D arrays have function of output size normalization. So convolution of 1D array can be solved like this:

Disadvantage of this solution is that the VI gets kinda slower since it have to work with 2D arrays.

Still I would appreciate some advice about how to find impulse response function with some more elegant way then fitting convolution.

[OpiCAN00]

I had some progress with Impulse Response Function PtByPt VI function. I calculate nice IRF when i split output array in the half and combine it vice versa. I was also able to find IRF with FFT after i round small values to the zero. But IRF from FFT is unstable for some combinations of sigma values and also Impulse Response Function PtByPt give unusable result,  when sigma of stimulus is higher than sigma of response.

I found some theory on that topic and I see it have to be investigate by Cross power spectra, but it is hard to follow without examples. https://www.sjsu.edu/people/burford.furman/docs/me120/FFT_tutorial_NI.pdf