I am trying to identify a system by means of its differential equation (i.e., Lapace representation). I put together a rather straightforward regression algorithm (similar to Proni's method for ARMA) under the assumption that the FFT of the system response is equivalent to the Laplace transform evaluated at \$ -j\omega \$ (with \$ \omega \$ restricted to \$ (-\pi , \pi ) \$, (but this doesn't seem to change anything).
It converges quite consistently but the poles I get are both positive and negative (i.e., stable and unstable if I interpret them as I intended in the s-domain) and both inside and outside the unit circle (again unstable and stable if I interpret them as a Z-transform instead).
I am no newby to DSP, but it is clear I have some misconception somewhere. Most likely on my interpretation of the FFT as simply a sampling of the Laplace transform on the imaginary axis or in the proper scaling of this axis. I have checked the equations back and forth multiple times, and I can't see what I am doing wrong.
The gist of it is simply to solve the equation:
$$ y(s)*(1 +a_1s + a_2s^2 + a_3s^3) = x(s)*(b_0 + b_1 s + b_2 s^2) $$
For the values of \$ a_i,b_i \$ that minimize the representation error, setting \$ s=j\omega\$ with \$ \omega \text{ in } (-\pi,\pi) \$ and making \$x(j\omega)\$ the FFT of a known input and \$ y(j\omega)\$ the FFT of the measured output (trimming the edge frequency bins as these are noisy).
I know that the FFT of a finite sampled sequence is simply sampling the Fourier transform and adding the aliased components into it (which in this case can be assumed to be zero, as the data stream is digitally filtered and subsampled).
I want to implement it in the continuous s-domain, as it is modeling an analog system and that would make the results easier to understand. So, before I throw all of this away and model it in the discrete Z-domain as an ARMA, can someone point to where my problem lies?