making a laplace s-domain plot from numerical data to decompose signal into decaying sinusoids

Question

I would like to make a 3D laplace s-domain plot from experimental data I have. The examples I have seen for this are when the function is already known and an analytical solution can be obtained (example How can I plot a 3D graph of a given Laplace Transform of a function?). I have experimental data (time series) where I don't know what functions are present, and am having difficulty generating the s-domain plot.

The approach I tried to take is

1. start with time domain signal y(t)
2. multiply it by e^(-sigma * t) for each value of sigma.
3. calculate the complex fourier transform
4. plot the s-domain, real vs imaginary. Look at poles and zeros.

I don't understand how to move from step 3 to step 4. The s-domain plot is a 3D plot with X, Y, and Z axes. I'm getting confused over what these different axes actually are.

X-axis = sigma, the value I am varying. Y-axis = the complex part of the output from step 3? Z-axis = ? magnitude of output from step 3?

Below is an example of what I'm trying to do in Python with an example signal that is similar to my experimental data. Ultimately I would like to do something similar to this question How to decompose a signal into exponentally decaying sinusoids?

import numpy as np
import matplotlib.pyplot as plt
from scipy import fftpack
#%% definitions

def synthetic_signal(t):
    #y = np.exp(-1.5 * t) * np.cos(40 * t) # simpler example
    y = 0.5 * np.exp(-10 * t) * np.sin(20 * 2 * np.pi * t + 0.1) + 1 * np.exp(-50 * t) * np.sin(250.0 * 2.0 * np.pi * t)
    return y

def laplace_function(function, t, sigma):
    y_exp = function * np.exp(-sigma * t)
    y_exp_fft = np.fft.rfft(y_exp)
    real = y_exp_fft.real
    imaginary = y_exp_fft.imag
    return real, imaginary
#%% look at result with a fixed sigma value to see if it makes sense.
t = np.linspace(0, 0.5, num=int(1e3), endpoint=False)
f = synthetic_signal(t)

l_real, l_imaginary = laplace_function(f, t, sigma=-1.5)

print('y_exp_fft =', y_exp_fft)

# plt.plot(l_imaginary, 'k.-')
# plt.plot(l_real, 'r.-')
plt.plot(l_real, l_imaginary, 'b.')

# then make a mesh plot
```

Answer 1

So I'm confused by your steps.

For me, the Laplace transform is usually calculated as $$ Y(s) = \int_{-\infty}^{+\infty} y(t) e^{-st} dt \tag{1} $$

So

start with time domain signal y(t)

that's OK, we have $y(t)$.

multiply it by e^(-sigma * t) for each value of sigma.

We can form

$$ y(t) e^{-st}$$

easily enough.

calculate the complex fourier transform

This is where I get confused: the formula (1) just does an integration, not a Fourier transform or FFT.

To do this step, I would just look at approximating the integral with a sum: $$ Y(s) \approx \sum_{t=0}^{T-1} y(t) e^{-st} $$ assuming the $dt$ is just 1.

The trouble with this approach is that I suspect you'll never get "poles" because the summation will always be finite (provided $y(t)$ is finite).

Answer 2

The steps given are correct given the relationship between the unilateral Laplace Transform and the Fourier Transform of a causal time function as demonstrated below:

Unilateral Laplace Transform:

$$X(s) = \int_0^\infty x(t) e^{-st}dt$$

Given $s= \sigma +j\omega$ as the real and imaginary components of the s-plane, this is:

$$X(s) = \int_0^\infty x(t) e^{-(\sigma+j\omega) t}dt= \int_0^\infty x(t) e^{-\sigma t}e^{-j\omega t}dt$$

Note that for $y(t)= x(t)e^{-\sigma t}$ That the unilateral Laplace Transform is just the Fourier Transform of $y(t)$ if $y(t)$ is causal:

$$Y(j\omega) = \int_0^\infty y(t) e^{-j\omega t}dt$$

Like the Fourier Transform, the resulting Laplace Transform is complex valued- every value returned for every possible value of $\sigma$ and $\omega$ has a magnitude and phase (when the Laplace Transform coverages). If we want to plot this two-dimensional surface, we can plot the magnitude and phase on two separate plots. In the linked posts by the OP, only the magnitude of the Laplace Transform was plotted. Given that, what the OP can do if it is desired to use the Fourier Transform for plotting is for each value of $\sigma$ multiply the time domain function by $e^{-\sigma t}$ and then plot the magnitude of the Fourier Transform (starting at $t=0$) for that vertical axis slice at $\sigma$ on the real axis. Note that the Laplace Transform only converges for $\sigma$ greater than the right most pole, so such a plot would start from the right most pole and scan over $\sigma$ values going to the right (increasing $\sigma$). This approach avoids having to compute the Laplace Transform directly but it will not reveal all the poles and zeros that are to the left of the rightmost pole (not in the region of convergence).

It is much easier to just make a plot from the resulting Laplace Transform in the case that it can be easily computed for a given time domain function. Given the Laplace Transform $X(s)$ as a function of s, simply compute the magnitude for each possible s on the 2D plane and plot that Surface (or phase if a phase plot is desired). This approach will show all the poles and zeros since the surface will be plotted regardless if the actual Laplace Transform converges there or not. In this case the surface viewed is that of the actual Laplace Transform only in the region of convergence.